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ABSTRACT 

A  material  model  accounting  for  strength  developed  earlier  for  two-phase  materials  is 
implemented  in  the  CTH  hydrocode.  The  strain  response  to  load  in  the  model  is  decoupled 
into  shear  and  volumetric  contributions  in  order  to  satisfy  the  model  implementation 
requirements  for  CTH.  Multi-phase  description  is  realised  via  constitutive  equations 
complementing  the  conservation  laws  for  a  material  represented  as  a  mixture  of  several 
phases.  Such  a  formulation  agrees  well  with  the  CTH  code  structure  and  is  suitable  for 
conventional  user  implementation.  The  implementation  has  been  applied  to  a  generic  material 
representing  sand  at  various  porosities.  The  constitutive  equations  and  equations  of  state  have 
been  fitted  in  order  to  describe  literature  data.  Numerical  illustrations  in  the  report 
demonstrate  agreement  of  the  calculation  results  with  the  anomalous  behaviour  observed  in 
the  literature  for  a  highly  porous  sand  at  shock  compression  and  a  good  description  of  the 
experiments  available  in  the  literature  on  the  explosion  of  a  sand-buried  charge. 
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CTH  Implementation  of  a  Two-Phase  Material 
Model  With  Strength:  Application  to  Porous 

Materials 


Executive  Summary 


Involvement  of  the  Australian  Defence  Forces  (ADF)  in  various  operational  theatres 
requires  improvement  of  protective  capability  against  buried  mines  and  Improvised 
Explosive  Devices  (IEDs).  The  threat  to  a  target  (e.g.,  a  vehicle)  due  to  a  detonating 
buried  device  is  twofold:  i)  the  structural  threat  due  to  mainly  gaseous  detonation 
products  that  are  loading  the  target  quite  slowly  in  the  millisecond  range  of  time  and 
over  a  wide  area  of  the  target;  and  ii)  the  impact  threat  due  to  the  soil  ejecta  colliding 
with  the  target  and  transferring  momentum  to  the  target  rather  quickly  within  the 
microsecond  range  of  time  and  localised  over  a  compact  area  of  the  target.  The  impact 
threat  is  the  most  immediate  and  important  one  before  the  structural  effects  take  place. 

To  develop  a  protection  capability  against  this  impact  threat,  an  enhanced  assessment 
of  the  response  of  porous  geological  materials  blanketing  a  detonating  device  and 
transferring  momentum  to  the  target  via  the  ejecta  impact  needs  to  be  undertaken. 
Such  an  assessment  can  be  performed  with  the  physics  modelling  tools,  hydrocodes, 
such  as  LS-DYNA  and  CTH  available  in  DSTO.  However,  even  with  the  use  of  these 
powerful  modelling  techniques,  behaviour  of  the  porous  geological  materials, 
commonly  sand  or  soil,  is  not  easy  to  describe  in  the  conditions  of  shock  loading. 
Complexity  in  the  behaviour  of  porous  materials,  demonstrated  for  example  in  [1],  is 
manifested  by  highly  non-linear  response  of  those  materials  due  to  their  multi-phase 
structure  with  drastically  different  compressibilities  of  constituents.  In  turn,  the  solid 
constituents  of  these  porous  mixtures  are  strength  resistant  and  strain  rate  sensitive. 
This  behaviour,  specifically  evaluation  of  the  parameters  responsible  for  the 
momentum  transfer,  is  not  well  predicted  within  the  traditional  approaches  using 
established  material  models  available  in  the  hydrocodes. 

In  order  to  improve  the  hydrocode  modelling  capability,  the  present  report  describes 
an  implementation  of  an  advanced  two-phase  material  model  [2]  that  takes  into 
account  the  multi-phase  nature  of  the  geological  materials  (sand,  soil,  etc).  Along  with 
the  air  contained  in  the  pores,  the  condensed  constituents  of  the  materials  are  also 
compressible  at  this  level  of  loads.  In  addition,  the  solid  constituents  are  strength  and 
strain  rate  sensitive.  In  the  present  report,  the  two-phase  model  [2]  is  implemented  in 
the  CTH  hydrocode  [3]  and  the  implementation  flowchart  is  briefly  outlined.  Sand  is 
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considered  as  a  model  material  in  the  test  calculations.  Mechanical  characteristics  of  the 
material  are  fitted  in  order  to  correlate  with  available  literature  data  for  the 
constituents,  namely,  the  air  and  quartz.  The  numerical  illustrations  demonstrate  a 
good  description  of  physical  features  that  are  typical  for  the  porous  materials  at  shock 
compression  and  are  difficult  to  describe  with  material  models  available  in  the 
hydrocode.  The  first  numerical  example  demonstrates  the  anomalous  behaviour 
observed  in  experiments  for  highly  porous  sand,  but  this  behaviour  is  not  customarily 
predicted  by  available  material  models  within  the  physics-modelling  framework  of  the 
hydrocode.  The  second  example  illustrates  formation  of  the  soil  ejecta  due  to  explosion 
of  a  buried  charge.  Comparison  of  the  numerical  results  with  available  literature 
experiments  shows  good  agreement.  The  illustrations  demonstrate  the  importance  of 
the  physics  modelling  for  the  description  of  parameters  responsible  for  the  momentum 
transfer  from  the  soil  ejecta  to  a  target. 

The  present  model  development  and  CTH  implementation  activity  is  also  performed 
as  a  part  of  the  joint  efforts  within  the  Modelling  and  Simulation  Focus  Area  of  the 
Conventional  Weapons  Technology  Group  (Terminal  Effects)  of  The  Technical 
Cooperation  Program  (TTCP). 
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1.  Introduction 


Neutralisation  of  the  effects  of  mines  and  buried  Improvised  Explosive  Devices  (IEDs)  on 
affected  vehicles  and  personnel  requires  an  analysis  of  the  target  loading  due  to  the  explosion 
of  a  buried  charge.  Such  an  analysis  involves  a  number  of  factors  to  consider.  The  most 
important  ones  amongst  them  are  the  size  and  area  distribution  of  the  charge,  mechanical  and 
thermal  properties  of  the  soil,  the  depth  of  burial,  and  the  target  stand-off  distance.  Other 
factors  such  as  method  of  initiation,  casing  arrangement,  and  the  properties  of  high  explosive 
are  still  important,  but  they  are  somewhat  secondary  for  charges  buried  under  a  layer  of  soil 
that  is  thick  enough,  as  soon  as  the  distribution  of  the  energy  to  be  released  at  a  given  depth  of 
burial  is  known. 

The  most  critical  factor  of  the  target  loading  analysis,  the  momentum  transfer  from  the  soil 
ejecta  to  a  target  such  as  a  vehicle's  floor,  is  mainly  controlled  by  the  density  and  velocity 
distributions  through  the  ejecta  thickness  when  the  ejecta  is  impacting  the  target.  The  ejecta 
impact  provides  the  highest  level  of  the  stress  transfer  contribution  into  the  target,  which  may 
possibly  be  transmitted  to  the  vehicle's  occupants  behind  the  floor  barrier  facing  the  ejecta. 
The  follow-up  blast  products  provide  a  comparably  smooth  and  gradual  loading,  highly 
energetic,  though,  to  affect  the  target's  structural  integrity  and  behaviour.  Another  important 
factor  for  the  analysis  is  separation  of  the  blast  products  from  the  soil  particles  (velocity  non¬ 
equilibrium  in  the  gas-solid  mixture).  This  factor  is  negligible  for  the  camouflet  blasts  (the 
detonation  products  are  contained  underground)  when  the  Depth  Of  Burial  (DOB)  is 
sufficiently  large  or  the  products  do  not  break  the  surface  before  the  ejecta  reaches  the  target. 
On  the  other  hand,  if  1)  DOB  is  not  sufficient  or/ and  2)  the  target  stand-off  distance  is  large, 
then  the  gaseous  products  may  break  through  and  the  gas-particle  separation  may  increase 
with  time  and  distance  of  travel  making  this  factor  more  influential.  At  the  same  time,  the 
localised  impact  damage  effects  on  the  target  due  to  the  charge  explosion  diminish  with 
increasing  stand-off  distance  or/ and  decreasing  DOB.  Moreover,  structural  effects  could 
possibly  be  analysed  while  neglecting  the  ejecta  impact  on  the  target  when  using,  as  a  basis 
for  the  analysis,  the  residual  energy  of  the  blast  products  after  having  broken  through  the  soil. 
Therefore,  the  scenarios  when  this  non-equilibrium  is  significant  are  of  less  interest  from  the 
viewpoint  of  the  localised  damage  effects  to  vehicle  and  vehicle  occupants.  When  evaluating 
these  effects,  the  mass-velocity  distribution  determining  the  momentum  transfer  is  believed  to 
be  the  more  relevant  ejecta-related  factor.  Thus,  omitting  the  structural  response  issue,  the 
present  consideration  is  restricted  to  the  velocity  equilibrium  case.  In  the  two-phase  analysis, 
synergistic  effects  might  also  be  taken  into  account  if  the  impact  and  blast  loadings  occur 
simultaneously. 

For  an  analysis  of  behaviour  of  porous  geological  materials  such  as  soil  or  sand  we  need  to 
pay  attention  to  a  quite  complex  material  response  requiring  an  advanced  approach.  The  first 
issue  of  the  analysis  is  the  porosity  of  the  substances.  The  second  issue  is  the  complex 
behaviour  of  the  solid  constituents  of  the  materials.  Specifically,  the  silicon  dioxide  (silica)  that 
is  the  most  common  component  of  the  solid  constituents  of  soils  is  known  for  by  its 
abundance  of  high-temperature  and  high-pressure  polymorphs.  Other  issues  involve  the 
presence  of  feldspar  minerals  and  water  (moisture  content)  that  may  easily  change  its  phase 
state  during  the  shock  loading.  As  a  starting  point,  we  consider  the  two-phase  representation 
of  sand  as  a  porous  material  with  the  solid  silica  constituent,  omitting  the  phase  transitions  of 
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silica  in  the  present  consideration.  Further  work  is  planned  to  account  for  a  phase  transition 
using  the  three-phase  modelling  approach  [1],  In  the  present  work  we  employ  a  two-phase 
material  model  with  strength  [4], 

Referring  to  the  model  implementation  process,  the  user-implementation  interface  of  modern 
hydrocodes  poses  certain  restrictions  on  the  implementation  of  material  models  in  codes  such 
as  LS-DYNA  [2]  or  CTFI  [3].  The  two  most  critical  restrictions  are  i)  the  use  of  the  fundamental 
form  of  the  conservation  laws,  namely,  the  mass,  momentum,  and  energy  conservation  laws, 
and  ii)  the  necessity  of  the  decoupled  representation  of  mechanical  response  in  the  form  of 
separate  spherical  and  deviatoric  stress  responses  of  a  material.  The  decoupled  response 
separates  the  pressure  contribution  calculations  for  the  momentum  update  within  the 
thermodynamic  block  of  a  code  from  the  deviatoric  stress  calculations  for  the  internal  material 
parameters  update  within  the  constitutive  block  of  the  code.  An  example  of  an 
implementation  dealing  with  these  issues  is  shown  in  [5]  for  a  rate  sensitive  strength  model. 
Having  specified  the  restrictions,  it  is  still  possible  for  code  developers  to  implement,  for 
example,  a  fully  coupled  system  for  the  full  deformation  gradient  or  stress  tensor  (e.g.  see  the 
manual  [3]  for  an  example  of  the  Transverse  Isotropic  (TI)  model).  The  code  developers  can 
also  escape  the  solitary  form  of  the  conservation  laws  and  implement  a  few  systems  of  the 
conservation  laws  for  material  components  (e.g.  see  [3]  for  a  few  implementations  of  the  Baer- 
Nunziato  multi-phase  flow  model).  However,  such  implementations  require  almost  full  access 
to  the  code  and  cannot  be  made  via  the  user-implementation  interface. 

The  initial  version  of  the  model  being  implemented  in  the  present  work  has  been  published  in 
[6].  Further  development  of  this  model  has  considered  the  inter-phase  heat  transfer,  which  is 
critical  for  porous  materials  and  allowed  the  author  to  simulate  the  anomalous  behaviour  at 
shock  compression  [7].  The  recent  version  [4]  of  the  two-phase  model  has  taken  into  account 
the  strength  resistance  of  condensed  phases  and,  in  turn,  rate  sensitivity  of  the  yield  limits  of 
the  phases.  The  present  work  deals  with  the  CTH  implementation  of  the  two-phase  model  [4] 
that  considers  the  inter-phase  heat  transfer  and  strength  effects  taking  rate  sensitivity  into 
account.  In  this  work,  the  CTH  implementation  is  realised  for  two-phase  porous  materials. 
Adjustment  of  the  implementation  to  any  two-phase  material,  where  phases  are  both 
condensed  and  the  phase  transition  is  possible,  is  a  routine  process  provided  that  the  phase 
transition  kinetic  is  known  (e.g.  such  an  adjustment  has  been  performed  in  an  in-house  wave 
code  [4]). 

When  implementing  the  model  in  a  hydrocode,  the  present  model  formalism,  similarly  to  [5], 
allows  one  to  decouple  stress  response  in  two  processing  streams:  1)  Elasto-Plastic  (EP)  block 
processing  Constitutive  Equations  (CE)  that  describe  the  evolution  of  deviatoric  elastic 
deformations  directly  linked  with  the  shear  stresses;  and  2)  Eulerian  Remap  Block  (ERB) 
processing  thermodynamic  relations  on  the  basis  of  Equations  Of  State  (EOSs)  that  describe 
the  volumetric  (bulk)  response  of  materials.  Strictly  speaking,  constitutive  equations  may  be 
used  for  calculation  of  not  only  the  deviatoric  stress  evolution,  but  other  internal  variables. 
Because  of  the  code  structure,  the  CE  block  is  used  at  the  code  point  when  deviatoric  stresses 
are  calculated  and  this  block  can  be  used  for  any  evolution  equations.  Therefore,  references  to 
the  CE  and  EP  blocks  are  interchangeable  in  this  report. 

The  procedure  of  implementation  in  CTH  is  briefly  described  below  with  references  to  the 
code  structure  available  in  the  public  literature.  The  main  CTH  structure  includes  the  Input 
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Block  (IB),  Lagrangian  Block  (LB),  and  Eulerian  Remap  Block  (ERB)  [8],  An  essential 
component  of  CTH  is  the  Database  Management  Block  (DMB)  that  arranges,  allocates,  stores, 
and  retrieves  the  data.  However,  this  block  is  not  separated  from  the  remaining  blocks, 
likewise  IB,  LB,  and  ERB  are  separated  from  each  other.  Rather,  DMB  code  fragments  precede 
every  call  of  any  major  subroutine  from  IB,  and,  particularly,  LB  and  ERB  blocks.  This  enables 
the  user  to  ignore  the  CTH  database  structure,  unless,  like  in  the  present  case,  the  database 
management  requires  additional  allocation  due  to  extra  requirements  for  internal  state 
variables.  Publications  [9,  10]  describe  the  code  structure  in  slightly  more  detail  with 
references  to  specific  routines  of  the  code.  The  present  implementation  includes  modifications 
to  i)  the  input  subroutines  of  the  EP  module  in  the  IB  part,  ii)  update  of  deviatoric  stresses  and 
extra  variables  in  the  LB  part;  and  iii)  update  of  key  thermodynamic  variables  and  extra 
variables  using  advected  parameters  and  EOS  in  the  ERB  part.  A  recent  CTH  implementation 
in  DSTO  [5]  describes  the  basic  steps  of  the  implementation.  However,  the  present  model  is 
essentially  more  complex  than  that  of  [5]  and  requires  more  extensive  use  of  extra  variables 
and  EOS  modules  of  CTH. 

To  illustrate  numerically  the  present  implementation,  a  response  of  generic  sand  is  modelled 
with  constituents  of  quartz  and  air.  The  model  uses  known  thermo-mechanical  data  for 
polycrystalline  or  fused  quartz  with  the  air  described  as  a  polytropic  gas.  Numerical  examples 
include  the  simulation  of  the  shock  wave  compression  of  a  highly  porous  sand,  and  the 
calculation  results  are  compared  with  the  shock  compression  data  for  silica  [11].  With  this 
calculation,  the  CTH  code  employing  the  present  implementation  demonstrates  its  ability  to 
describe  the  anomalous  behaviour  of  the  Hugoniots  [11]  at  a  high  porosity.  Another  example 
is  the  simulation  of  the  expansion  of  an  initiated  High  Explosive  (HE)  charge  buried  under  a 
layer  of  sand.  The  experiments  [12]  are  compared  with  the  CTH  calculations  using  two  CTH 
database  models  and  the  implemented  model.  The  numerical  results  show  that  the  present 
model  demonstrates  an  improvement  in  predictions  when  comparing  with  the  models 
available  in  CTH. 


2.  Constitutive  Model 

2.1  Conservation  Laws  and  Constitutive  Equations 

The  system  of  equations  of  the  two-phase  model  with  strength  contains  three  conservation 
laws  for  mass,  momentum  and  total  energy  with  the  stress  tensor  in  the  decoupled  form  a,y  =  s,y 
-  p5ij  (see  [4]): 


dp  |  dpuk  _Q 
dt  dxk 

dp[e  +  u  / 
dt 


dpUj  |  d(puiuj-sj)  |  dp  _Q 
dt  dx ,  dx, 


Puj\e  +  u/2)-sijui+  puj  j 


(1) 


dx, 


=  0. 


Denotations  of  the  physical  variables  are  standard  and  described  in  [4],  The  strength  response 
is  described  by  'specific'  strain  e,j  when  small  shear  deformations  are  assumed  (see  [1, 4,  5]). 
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Description  of  the  strength  response  uses  the  same  symmetrical  velocity  gradient  as  the 
majority  of  other  strength  sensitive  models  in  CTH.  The  response  is  described  by  the 
following  system  of  equations  [4]: 


dpe„  dpe,uk  d 


c  ■>/ 


dUk/ 


dt 


dx. 


dx . 


dx: 


■  +  ■ 


dx,. 


•  s,j  =  -  +  (i  -  #k!2) )  -  p *  jv  > 


(2) 


where  yyrl)  =  p  [e,;  +  (1  -  0)  A y]/ n  and  cpp1)  =  p  [ty  -  0  A,,]/ n.  Superscript  indices  refer  to  the 
material  constituent  number  and  n  and  ti  are  the  relaxation  time  functions  obtained  for  rate 
sensitive  materials  of  the  phases  from  the  corresponding  yield  limits  versus  strain  rate  (see 
[13]).  The  parameters  c  and  0  are  mass  and  volumetric  concentrations  of  the  first  phase  and  Ay- 
are  components  of  the  strain  inter-phase  imbalance  [4],  which  evolves  in  accordance  with  the 
following  constitutive  equation  [4]: 


dp^j  ^  dpXijuk 
dt  dxk 


)• 


(3) 


Function  ip  in  the  last  term  of  the  right-hand  side  of  the  equations  (2)  is  the  compression 
kinetic  from  the  following  constitutive  equations  describing  the  phase  transition  and 
compaction  due  to  relative  compressibility  between  the  phases  [4]: 


dpc  dpcuj 
dt  dx  j 


pep  =  m0  , 


dpQ  dpOuj 
dt  dx  j 


PV  ■ 


(4) 


The  inter-phase  heat  transfer  is  described  with  the  following  equation  for  the  entropy 
disequilibrium  j  [4]: 


dp%  ,  dpx u 
dt  dx 


pco  . 


(5) 


The  inter-phase  heat  exchange  kinetic  co  is  also  to  be  specified  [4,  7], 

The  system  of  differential  equations  of  this  subsection  describes  the  behaviour  of  a  two-phase 
mixture.  However,  this  system  calculates  only  the  evolution  of  independent  thermodynamic 
and  kinematic  variables  (with  the  exception  of  entropy  described  indirectly  by  the  energy 
conservation  law).  Therefore,  the  system  should  be  completed  with  algebraic  equations  for 
dependent  thermodynamic  variables,  which  includes  the  kinetic  functions  for  the  constitutive 
equations  and  an  equation  of  state. 


2.2  Gibbs  Energies  and  Thermodynamic  Potentials  of  two-phase  mixture 

In  order  to  specify  relations  between  the  dependent  variables  e,  p,  Sy„  and  T  and  independent 
variables  p,  c,  0,  S,  cy.  Ay,  and  y,  we  need  to  recall  the  relations  between  the  averaged  variables 
and  the  variables  associated  with  individual  phases  [4].  The  phase-specific  variables  are 
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marked  by  superscripts  ‘V  and  '2'  in  parentheses  for  the  first  and  second  phases,  respectively, 
with  the  first  phase  being  air  in  the  case  of  a  porous  material.  Specifically,  the  phase 
independent  variables  are  calculated  from  the  averaged  variables,  using  the  concentration  and 
disequilibrium  parameters  from  the  following  relations  [4] : 


p(i)  =  P  c/9 ,  pi 2)  =  P(1  -  c)  /(I  -  0 ) ,  SO)  =  Vi(S  +x)/c,  S( 2)  =  1/2  (S  -x)  /(I  -  c) , 


df)  =  9  (e,j  +  hj  (1  -  6 ))  /c ,  ei/2)  =  (1  -  6 )  (e,,-  -  Ay  0)  /(I  -  c) . 


The  dependent  averaged  variables  are  calculated  via  the  phase  ones  as  follows 
e  =  c  eO)  +  (1  -  c)  e(2),  p  =  p2eP  =  9  pi1)  +  (1-9)  pi2),  T  =  es=(  TO)  +  T(2))/2  . 
Sy=eej'=(ce"+(l-c)e*%=eSv+(l-e)Se\ 


(6) 


(7) 


The  fundamental  thermodynamic  characteristics  used  for  the  description  of  the  mixture 
equilibrium  is  the  Gibbs  energy  for  each  phase  pu  =  eik)  +  pik)/ pik)  -  Tik)Sik )  -  c1?O0s1;D,  fc  =  1,2. 
Affinity  of  the  Gibbs  energies  characterizes  phase  stability  (see  [1]  for  details  and  an  example). 
Similar  generalization  of  the  Gibbs  energy  for  specific  materials  has  been  considered  earlier 
for  two-phase  mixtures  without  strength  in  [6, 14]  and  for  the  mixtures  with  strength  in  [1, 4, 
5].  Using  the  Gibbs  energy,  the  thermodynamic  potentials  of  the  mixture  are  (see  [4]  for 
details) 


„0>  „(2) 

A  =  =  e<»  +  ^  - r»>S">  -  4»4'>  -  +  T-e’S®  +  =  ft  -  ft2 , 

Q„=eli(=e( 


n  =  =  -  - 


(1)  + 


-  +  ■ 


,(1)e»y))(J?  ~sij)~P{1)  -(p(2)42)X42) 


~  ^  ijsij  =n  o-AijS  = 


(8) 


T  =  ez  =(t(1>-T(2>)/2. 


P 


In  order  to  satisfy  thermodynamic  correctness,  the  constitutive  relations  have  the  following 
restrictions  on  the  right-hand  sides  (see  [4]): 

cp  =  A-cpo  ,  ip  =  Tlo-ipo  ,  co  =  W-qo  -  (9) 

Thus,  when  equations  of  state  are  given  for  each  phase  in  the  form 

e(i)  =  e(i)(p(i),  eif1),  si1))  ,  ei2)  =  ei2)(pi2),  eif2),  si2))  ,  (10) 

the  conservation  laws  (1)  and  the  constitutive  equations  (2)-(5)  generate  a  closed  system  of 
equations.  The  system  can  be  solved  numerically  if  functions  ei1)  ,  ei2)  and  <po,  </)(>,  coo  are 
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specified.  An  example  of  such  a  specification  for  porous  sand  will  be  shown  in  the  next 
Section. 


3.  Constitutive  Relations  and  EOS 


3.1  Equation  of  State  of  two-phase  mixture 

For  porous  materials,  we  specify  the  first  phase  to  be  air  and  the  second  phase  a  solid 
constituent.  In  the  present  implementation,  to  determine  EOS  (10)  for  the  condensed 
constituent  (the  second  phase  in  a  porous  mixture  denoted  by  superscript  '2'),  we  choose  a 
reduced  form  of  the  Mie-Gruneisen-type  EOS  developed  in  [15]  and  used  in  [1,  4-  7]: 


2a  „ 


■(S‘-  -if  +2plblSf-*2d  +cvsT„S'- 


+  Zc_Zo 
A  P 


(li) 


Here,  ao,  fio,  and  pi  (Gruneisen  coefficient)  are  material  constants,  cvs  is  the  thermal  capacity, 
and  5  =  p/pos  with  initial  density  pos-  The  constant  ao  is  the  bulk  sound  velocity  that  is  linked 
with  the  longitudinal  and  shear  sound  velocities  do  and  bo  as  follows 

ao2  =  do2  -  4b02/3  ,  (12) 


and  d  in  (11)  is  the  second  invariant  of  the  'specific'  strain  deviator:  d  =  e,j  •  ep  / 2.  It  should  be 
noted  that  the  arguments  and  the  function  itself  for  the  internal  energy  (11)  refer  only  to  the 
second  constituent  of  the  two-phase  mixture. 

Shear  stresses  are  calculated  from  (11)  for  the  solid  phase  as 


s-  =  e„ 


=  e,,d  =  e,e.. 

d  e,  d  ij 


=  2  p{ 


\b\S>"\=lpblSf'pet,  =2  Gpe, 


(13) 


Where  the  shear  modulus  G:  G  =  pb^S11" . 


For  the  first  constituent  in  porous  material,  the  air,  EOS  in  (10)  is  selected  in  the  ideal  gas 
form: 


r(r  + 1) 


(d’-l ,)+c,V' 

exp 

fs/  ) 

-1 

K/  ) 

(14) 


Here  y  is  the  polytropic  gas  exponent,  Co  is  the  sound  velocity,  cvg  is  the  thermal  capacity  of  the 
gas,  and  5  =  p/pog  with  the  initial  density  pog-  The  following  natural  restriction  follows  from 
the  definition  of  the  ideal  gas  EOS:  Co2  =  y  (y  +  1 )  cvgTo.  In  this  case,  the  variables  relate  only  to 
the  phase  '1'.  The  initial  densities  of  the  two  constituents  pog  and  pos  are  used  along  with  the 
initial  density  of  porous  material  po  in  order  to  define  initial  mass  and  volumetric 
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concentrations  c  and  6  from  (6).  These  concentrations  and  density  are  used  as  initial  data  for 
the  system  of  conservation  laws  and  evolution  equations. 


3.2  Constitutive  Relations  characterising  strength 

The  accounting  for  the  strength  of  the  second  constituent  (in  the  case  of  porous  material)  is 
realised  via  the  constitutive  equations  (2).  In  order  to  close  the  equations,  the  relaxation  time 
functions  for  solid  constituents  with  strength  should  be  chosen.  In  the  present  case,  the 
function  for  the  second  constituent  is  chosen  in  the  following  form  [5, 13]: 


t(p,s,T)  =  t0 


exp 


(D0  +  He), 


N0  +  Me 


(15) 


where  s  determined  from  s2  =  s,j  -s ij  is  the  second  invariant  of  the  stress  deviator.  In  the  present 
work  as  well  as  in  [13],  the  variable  e  is  associated  with  the  elastic  portion  of  deformation  e  = 
s/  (2 G)  and  to.  Do,  H,  No,  M  are  material  constants.  It  should  be  noted  that  the  constants  in  (11), 
(14)  and  (15)  are  specific  to  each  phase  and  in  the  present  implementation  the  function  (15)  is 
relevant  only  to  the  second  phase  of  porous  mixture.  Thus,  using  the  function  (15)  for  the 
second  phase  as  T2,  we  can  determine  cpifi2)  and,  assuming  absence  of  any  strength  effects  in  the 
first  gaseous  phase,  we  take  cp,/1)  =  0  in  (2)  and  (3)  for  the  present  case. 


In  order  to  facilitate  choice  of  constants  in  (15),  an  algorithm  of  fitting  the  constants  suggested 
in  [13]  has  been  realised  in  the  present  implementation.  Let  us  briefly  outline  the  algorithm. 
Firstly,  the  user  has  to  select  two  yield  limit  points  Yi  and  Y2  at  different  strain  rates  ex  and 
e2 .  It  should  be  noted  that  this  choice  is  a  responsibility  of  the  user  and  the  present  fit  does 
not  guarantee  a  proximity  of  any  pre-determined  experimental  or  hypothetical  curve  passing 
through  the  two  given  points  in  the  (Y-  e )  plane  to  the  curve  obtained  with  this  fit  except  for 
the  two  points  specified.  Secondly,  the  user  needs  to  select  constants  No  and  M  as  an  initial 
density  and  a  multiplication  factor  of  defects  in  the  condensed  material.  As  shown  in  [16], 
these  constants  have  a  much  smaller  influence  on  the  fit  of  the  curve  than  the  constants  To  and 
Do.  The  parameter  H  is  tabulated  for  future  use  to  describe  the  hardening  behaviour  of 
materials  and  cannot  be  used  for  modelling  such  a  behaviour  with  the  present  choice  of  £  in 
(15)  (see  discussion  on  this  topic  in  [5]).  Therefore,  it  is  taken  H  =  0  in  the  present  case.  The 
idea  of  fitting  constants  To  and  Do  in  [13]  is  based  on  an  approximation  of  the  stationary 
solution  of  the  viscoelastic  model  equations  [17],  It  was  observed  in  [13]  that  this  solution  is 
sufficiently  close  to  the  yield  limit  point.  The  corresponding  stress  relaxation  equations  for 
characterising  the  plasticity  state  (a  single-phase  analogue  of  the  equation  (2))  are  as  follows 
(see  [5]): 


dt 


du,  du  , 


-  +  ' 
dx ,  dx; 


£'e 


(16) 
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In  the  uniaxial  stress  conditions  (022  =  0),  the  yield  limit  is  Y  =  On  and  s  =  dul  /ax,  .  Therefore, 
the  stationary  point  of  (16)  after  index  summation  can  be  found  from 


2Gs  =  — . 

T 


The  last  equation  with  r  substituted  by  (15)  results  in 
.  .  .  Y  .  Ta  Dn 


2  G 


N  ' 


Y 


[3  . 

f  MY  ) 

—  +  In 

v  2 

l  2  Ng  2  Gj 

(17) 


that  allows  one  to  find  To  and  Do  when  two  pairs  {Yi,  a, }  and  {Y2,  s2 }  are  given.  Usually,  these 
data  are  known  from  the  Split  Hopkinson  Pressure  Bar  (SHPB)  tests  at  compression. 
However,  the  same  algorithm  can  be  used  for  the  tensile  tests  as  well.  For  simplicity,  we 
assume  that  when  necessary  we  can  utilize  another  pair  at  tension  with  yield  limits  Yu  and  Y> t 
using  the  same  strain  rates.  In  the  present  implementation  we  use  p  <  0  as  a  condition  for  use 
of  the  tensile  yield  limits. 


3.3  Thermo-Mechanical  and  Mass-Exchange  constitutive  relations 

Closure  of  the  constitutive  equation  (5)  describing  inter-phase  heat  exchange  can  be  done  via 
specification  of  the  right-hand  side  co  for  the  equation.  In  conjunction  with  the 
thermodynamics  condition  (9),  the  required  inter-phase  heat  exchange  kinetic  coo  in  (5),  (9)  is 
used  in  the  following  form  [4] 


COg  =  2(0  g 


h  Til)+T(2) 
p  T(1]T(2) 


h  =  To 


neff 

dl 


7  _  _ ^1^2 _ 

eff  ~  k2o  +  kfy-o)  ’ 


(18) 


where  h  is  the  heat  transfer  coefficient,  do  is  a  typical  material  structure  size  (e.g.  the  grain  or 
cell  size  for  porous  materials).  As  is  a  dimensionless  coefficient  (a  function  of  the  factors 
promoting  or  resisting  the  heat  exchange,  see  [7]  for  details;  in  the  present  work  As  =  400),  and 
ki  and  ki  are  the  thermal  conductivities  of  the  phases.  Initial  entropy  concentration  j  can  be 
calculated  from  initial  temperatures  TW  =  dc(r> / d s(1  \  T(2)  =  3e(2) / ds,2  ]  and  from  densities  (6)  if  the 
temperatures  initially  are  not  in  equilibrium.  Omitting  this  rather  exotic  case  in  the  present 
implementation,  we  take  s  =  0  and  j  =  0  as  the  initial  data  (the  initial  temperatures  of  the 
phases  take  the  reference  values  in  this  case). 

For  closure  of  the  mass  exchange  constitutive  equation  (the  first  equation  from  the  pair  (4)), 
the  function  cp  should  be  defined  via  the  specification  of  cpo  from  (9).  In  the  case  of  a  phase 
transition  occurring  between  two  phases  of  a  two-phase  material,  an  example  of  the  definition 
of  this  function  can  be  found  in  [4].  In  the  present  case  of  a  porous  material,  absence  of  the 
inter-phase  mass  exchange  can  be  assumed.  Therefore,  the  choice  cp  =  0  in  (20)  for  the  porous 
mixture  finalises  definition  of  the  mass-exchange  kinetic. 
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3.4  Compaction  kinetic 


The  last  equation  to  be  specified  is  the  right-hand  side  ip  for  the  constitutive  equation 
describing  the  inter-phase  compressibility  (the  second  equation  from  pair  (4)).  This  function 
(the  compaction  kinetic)  is  normally  determined  from  the  SHPB  compaction  tests  when  the 
porous  sample  is  confined  within  an  assembly  maintaining  a  constant  pressure  (e.g.  see  [18]) 
or  within  a  container  with  rigid  walls  maintaining  zero  lateral  strain.  In  the  assumption  of 
uniaxial  strain,  one-dimensional  homogeneous  deformation  (u  =  t  -x  +uo),  and  absence  of 
mass  exchange,  the  system  of  equations  for  determination  of  the  compaction  kinetic  is 
reduced  from  (1-5)  to  the  following 


dlnP  |  a_q  den  2  {l-OW 

dt  dt  3  p  p 


de21 

dt 


P 


dA  22 
dt 


d0_ 

dt 


-v. 


A  uY  ■> 


dA  j  j 
dt 


P 


(19) 


Here,  £  is  the  rate  of  deformation,  which  can  be  a  function  of  time,  and  the  entropy 
dissipation  term  in  (19)  is  (see  [4]): 


<D  =  — 
T 


(pA  +  y/U0  +a>x¥  + 


6<p. 


(i)  „(i) 
ij  *ij 


+ 


(i  -eWs? 


p 


(20) 


In  the  absence  of  inter-phase  mass  exchange  for  the  present  case,  cp  =  0  and,  because  of  the 
gaseous  first  phase,  cpd1)  =  0.  Keeping  in  mind  that  the  relaxation  function  ti  is  assumed  to  be 
determined  from  appropriate  SHPB  data  for  the  second  constituent,  we  can  rewrite  (19)  in  the 
following  further  reduced  form 


d\np 

dt 


=  0, 


d6_ 

dt 


-V  > 


-2  G2s 
3 


dS_ 

dt 


=  0 


r2  ’  dt 


--2  G2£ 
3 


f 

y/  no  +  tfAF  + 

v 


(1  -6)<P? 


P 


‘J  M) 
^  ij 


(21) 


where  Sif2)  =2Gi(eij  -  X,y  6).  Therefore,  pA,/  =  -  s,/2)/2G2,  and  pe,j  =  (1-0)  The  system  (21) 

for  the  stress  tensor  has  been  obtained  from  (19),  when  multiplying  the  corresponding 
equations  by  G2  (assuming  this  modulus  to  be  a  constant  within  a  single  time  step)  and  by 
manipulation  of  the  'specific'  strain  and  the  strain  disbalance  to  produce  s,/2).  Thus,  the  system 
(21)  allows  one  to  calculate  p,  0,  j,  S,  en,  £22,  An,  and  A22  as  soon  as  the  right  hand  sides  in  (21) 
are  given.  Here  c  is  constant  due  to  the  absence  of  mass  exchange  in  the  porous  mixture.  The 
stresses  a,,  =  s /,  -  p5,,  are  computed  from  (7)  as  follows 

ail  =  (1-0)  S1T2)  -  0  pd)  +  (1  -  0  )  p( 2)  ,  U22  =  (1-0)  S22 (2)  -  0  p(!)  +  (1  -  0 )  p(2)  ,  (22) 
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where  pCO  and  p(2)  are  obtained  from  EOSs  (14)  and  (11),  using  the  thermodynamic  identity  for 
individual  phases,  similarly  to  (7).  Summarizing,  while  the  functions  ti  and  co  having  been 
pre-determined,  the  function  ip  can  be  fitted  from  the  stress-strain  curves  recorded  during 
SHPB  compaction  experiments. 

As  a  first  approximation  case  for  the  compaction  kinetic  of  sand,  the  function  ipo  in  (9)  was 
taken  from  [6]  in  the  following  form: 

ipo  =  ipood  (1-0)  (0  -0o(p)) ,  at  0  >  0o(p)  , 

(23) 

ipo  =  0  ,  otherwise, 

where  0o(p)  is  arranged  in  such  a  way  that  0o  =  0oo  at  p<pc  and  0o  quickly  drops  from  the  initial 
value  0oo  to  zero  in  a  smooth  way  when  p  exceeds  pc  (a  simple  polynomial  approximation  with 
the  power  exponent  nc).  This  function  (23)  is  used  in  the  present  calculations  as  a  compaction 
kinetic. 


4.  CTH  Implementation 

The  present  implementation  affects  the  CTH  blocks,  IB,  LB,  and  ERB,  mentioned  earlier,  and 
allocates  extra  memory  in  a  DMB  module  of  LB  for  the  Internal  State  Variables  (ISVs). 


4.1  Input  Block 

For  the  input  block,  IB,  three  subroutines  [5,  9,  10,  19],  UINEP.FOR,  UINISV.FOR,  and 
UINCHK.FOR,  are  modified.  In  addition,  modifications  of  two  substitute  subroutines, 
EOSVEI  and  EOSVEK,  of  the  EOS  input  module  are  introduced  in  two  subsequent  subsections 
as  described  below. 

4.1.1  EOSVEI  modifications 

In  order  to  conduct  the  EOS  calculations,  the  substitute  subroutines  available  for  an  EOS  one- 
component  VE  model  are  replaced  by  subroutines  representing  the  present  EOS  combined 
from  (11)  and  (14)  via  the  mixture  rule  (7).  Thus,  the  'VE'  identifier  along  with  the 
corresponding  input  data  refers  to  the  present  EOS  model.  The  input  EOS  data  are  taken  from 
the  EOS_data  array  by  selection  of  the  material's  name  in  the  corresponding  section  related  to 
the  specified  model.  However,  this  dataset  has  to  be  replaced  by  explicitly  given  data  that  are 
necessary  for  the  present  model.  In  doing  so,  any  material  name  can  be  used  as  the  dummy 
name  in  the  EOS  input  (for  example,  the  material  name  'pmma'  is  used  in  the  input  decks 
below). 

The  EOS  subroutines  require  the  initial  densities  of  the  mixture  and  phases  (they  are 
calculated  later  on  in  the  EP  part  of  the  input).  Thus,  the  initial  density  of  the  porous  material 
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and  the  densities  of  the  first  and  second  phases  have  to  be  calculated  beforehand  and 
introduced  explicitly  via  three  successive  variables  that  are  an  input  in  the  original  VE 
subroutine  (they  are  denoted  by  'RO',  'TO'  and  'CV',  respectively,  by  the  subroutine 
convention). 

4.1.2  EOSVEK  modifications 

The  extra  variables  to  be  used  in  the  EOS  block  are  introduced  in  this  subroutine  because  the 
EOS  substitute  subroutines  do  not  provide  access  to  the  extra  variables  defined  in  the  CE 
block.  Thus,  the  value  exchange  between  the  EOS  and  CE  extra  variables  at  a  common  point  of 
the  code  is  needed  (see  subsection  4.3.1  below).  The  extra  variables  required  by  the  EOS 
thermodynamics  block  for  use  with  the  model  are:  i)  symmetric  'specific'  strain  5-set  tensor 
(ER  variables);  ii)  scalar  mass  concentration  c  (CRMS  variable);  iii)  scalar  volume 
concentration  0  (RTET  variable);  iv)  scalar  entropy  disequilibrium  j  (XIR  variable);  and  v) 
symmetric  strain  disbalance  5-set  tensor  Ay  (LMR  variables).  Specification  of  the  internal 
parameters  for  each  of  the  variables  is  described  in  more  detail  below,  when  outlining 
modifications  for  the  subroutine  UINISV.  Along  with  the  thermodynamic  variables  provided 
by  CTH,  this  set  is  sufficient  for  calculation  of  all  necessary  thermodynamic  functions  in  the 
EOS  module. 


4.1.3  UINEP  modifications 


UINEP.FOR  modifications  read  in  the  data  from  VP_data  input  file  into  the  VPUINP  array 
allocated  for  the  EP  related  input  data  [5, 19] .  Because  of  the  access  restrictions  to  the  EOS 
user-implementation  block,  both,  EOS  and  CE  data,  are  taken  from  the  VP_data  input  file  that 
results  in  a  fairly  significant  block  of  data.  The  data  needed  for  the  model  input  are:  'RHO', 


'COG',  'GAMG',  'RHOS',  'COS',  'BOS',  'ALF',  'BET',  'GAMS',  'CVS'/PRSO',  'LGEP1',  'LGEP2', 
'YIC',  Y2C',  'Y1T',  Y2T',  'AN2', ' AM2',  'TETOO',  'PCR',  'ANCR',  'PSIOO',  'DO', ' AKG',  'AKS', 
and  'AHTS'. 


These  constants  represent  the  following: 

initial  density  of  porous  material  for  RHO  (po); 

the  bulk  sound  velocity  of  the  gaseous  phase  for  COG  (co  in  (14)); 

the  polytropic  gas  exponent  for  GAMG  (p  in  (14)); 

initial  density  of  the  solid  phase  for  RHOS  (pos  in  (11)) 

the  bulk  sound  velocity  of  the  solid  phase  for  COS  (flo  in  (11)); 

the  shear  sound  velocity  of  the  solid  phase  for  BOS  (bo  in  (11)); 

the  bulk  modulus  exponent  of  the  solid  phase  for  ALF  (flo  in  (11)); 

the  shear  modulus  exponent  of  the  solid  phase  for  BET  (J3 o  in  (11)); 

the  Gruneisen  coefficient  of  the  solid  phase  for  GAMS  (po  in  (11)); 

the  thermal  capacity  of  the  solid  phase  for  CVS  ( cvs  in  (11)); 

ambient  pressure  for  PRSO  (po  in  (11); 

decimal  logarithm  of  strain  rate  related  to  the  first  test  point  for  LGEP1  (log  sx  in  (18)); 
decimal  logarithm  of  strain  rate  related  to  the  second  point  for  LGEP2  (log  A,  in  (18)); 
the  yield  limit  at  compression  related  to  the  first  point  for  YIC  (Yi  in  (18)); 
the  yield  limit  at  compression  related  to  the  second  point  for  Y2C  (Yi  in  (18)); 
the  yield  limit  at  tension  related  to  the  first  point  for  Y1T  (Yit  in  (18)); 
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the  yield  limit  at  tension  related  to  the  second  point  for  Y2T  (Yit  in  (18)); 

the  initial  defect  density  related  to  the  solid  phase  for  AN2  (No  in(15)); 

the  defect  multiplication  coefficient  related  to  the  solid  phase  for  AM2  (M  in(15)); 

the  starting  volume  concentration  for  TETOO  (0oo  in  (23)); 

the  critical  pressure  of  the  compaction  kinetic  for  PCR  (pc  in  (23)); 

the  power  exponent  of  polynomial  approximation  of  the  kinetic  for  ANCR  (nc  in  (23)); 

the  proportionality  constant  of  the  compaction  kinetic  for  PSIOO  ( ipoo  in  (23)); 

the  characteristic  material  structure  size  in  the  heat  transfer  kinetic  for  DO  (do  in  (18)); 

the  thermal  conductivity  of  the  gaseous  phase  for  AKG  (ki  in(18)); 

the  thermal  conductivity  of  the  solid  phase  for  AKS  (fe  in(18));  and 

the  dimensionless  coefficient  of  the  heat  exchange  kinetic  for  AHTS  (As  in  (18)). 

It  should  be  noted  that  the  model  uses  non-standard  input  units  in  cm  (length),  g  (mass), 
lOpsec  =  10_5sec  (time),  and  °K  (temperature).  The  strain  rate  £  for  LGEP1  and  LGEP2  is 
taken  in  inverse  seconds.  If  the  value  Y1T  is  negative  then  the  data  Y1T  and  Y2T  are  taken  to 
be  identical  to  the  set  Y1C  and  Y2C.  The  derived  pressure  unit  in  this  dataset  is  GPa.  At  the 
end  of  the  modifications,  the  initial  value  of  Poisson's  ratio,  and  the  bulk  and  shear  modulus 
are  calculated  and  initially  checked  in  the  same  UINEP.FOR  subroutine. 

4.1.4  UINCHK  modifications 

UINCHK.FOR  modifications  call  new  subroutines  VERCHK.FOR  and  VERCHX.FOR.  Before 
the  calls,  the  modified  code  checks  if  the  assigned  MODLEP  number  for  the  present 
constitutive  model  is  in  agreement  with  the  substituted  EOS  number  MEQ  [19].  The  next  part 
of  the  modifications  assigns  the  necessary  model-type  attributes  of  the  model  such  that 
deviatoric  stresses  will  be  calculated.  Then,  a  standard  call  to  the  subroutine  SI2CTH 
introduces  the  unit  transformation  constants  into  a  part  of  the  input  array  VPUINP.  The  new 
subroutine  VERCHK.FOR  transforms  the  constants  into  the  CTH  units  from  the  non-standard 
input  units,  introduces  global  constants  (the  'GC'  part  [19]  of  the  array  VPUINP)  such  as  the 
ambient  pressure  and  reference  temperature  used  in  (11)  and  fills  in  the  'DC'  part  of  the  array 
with  several  derived  constants  used  for  EOS  and  CE  calculations  such  as  the  time  relaxation 
constants  for  (15). 

The  new  subroutine  VERCHX.FOR  called  afterwards  replaces  the  EOS  input  data  (UI)  by  the 
data  from  the  VP_data  file,  and  recalculates  auxiliary  constants  for  the  initial  cycle  of  the  EOS 
calculations  (before  the  EP  input)  filling  in  the  tail  of  the  UI  array  for  the  EOS  block  [19]  with 
these  constants.  The  EOS  constants  previously  taken  from  the  EOS  analogue  of  the  VP_data 
file  are  now  replaced  by  the  data  described  above  in  the  present  section,  which  is  done  when 
calling  VERCHX.FOR.  The  whole  set  of  variables,  when  comparing  with  only  3  variables 
introduced  in  EOSVEI.FOR,  is  prepared  for  the  EOS  module  in  this  subroutine.  However,  the 
mass  calculation  is  processed  earlier  than  the  elasto-plastic  input  and  should  be  preceded  by 
input  of  the  EOS  data.  As  mentioned  in  subsection  4.1.1,  the  code  uses  the  EOS  data  for  the 
initial  cycle  calculation  of  mass  arrays.  In  order  to  use  the  input  data  correctly,  when  using  the 
EOS  identifier  in  the  EOS  section  of  the  CTH  input,  we  explicitly  assign  values  of  the  initial 
mixture  and  phase  densities  for  the  EOS  model  to  be  substituted,  while  leaving  the  remaining 
parameters  of  the  EOS  input  intact.  These  values  should  be  in  agreement  with  the  value  of 
'RHO'  in  the  VP_data  input,  and  with  the  values  for  the  density  of  the  gaseous  phase 
calculated  from  COG,  GAMG  and  PRSO  and  for  density  RHOS  for  the  solid  constituent. 
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4.1.5  UINISV  modifications 

UINISV.FOR  modifications  make  two  calls  for  a  standard  subroutine  MIGSEX  [19]  setting  up 
default  values  for  extra  variables  and  for  a  new  subroutine  VEREXV.  The  subroutine  VEREXV 
sets  up  extra  variables  to  be  used  in  the  CE  block.  Some  values  of  the  extra  variables  are 
shared  via  the  VERSWP  subroutine  outlined  below  in  the  ERB  modifications.  This  value 
exchange  is  needed  in  the  present  implementation  because  access  to  the  EOS  block  is  limited 
and  extra  variables  defined  in  EOSVEK  above  cannot  be  easily  reached.  The  whole  set  of  extra 
variables,  which  is  required  by  both  the  CE  and  the  EOS  thermodynamics  blocks  for  use  with 
the  model  and  being  specified  in  the  VEREXV  subroutine,  is:  i)  symmetric  'specific'  strain  5- 
set  tensor  e,y  (E  variables);  ii)  scalar  mass  concentration  c  (CMS  variable);  iii)  scalar  volume 
concentration  0  (TET  variable);  iv)  scalar  entropy  disequilibrium  j  (XI  variable);  v)  symmetric 
strain  disbalance  5-set  tensor  A y  (LM  variables);  vi)  the  old  density  at  the  start  of  the  time  cycle 
calculation  (RHOO  variable);  and  vi)  entropy  at  the  start  of  the  time  cycle  calculation  (ENT 
variable).  This  set  is  sufficient  to  calculate  all  the  necessary  thermodynamic  functions  in  both 
the  EP  and  EOS  modules.  The  tensor  character  for  variables  E  and  LM  and  the  scalar  character 
for  the  remaining  variables  are  specified  by  proper  selection  of  the  parameter  ITYPE  in  the 
subroutine.  All  the  variables  are  being  advected  at  the  Eulerian  remap  step  except  for  the  last 
two,  which  are  demanded  by  appropriate  selection  of  the  parameter  IADVCT  for  the  variables 
defined  above  (see  [19]).  The  dimensioning  parameter  (array  RDIM  [19])  is  set  up  in 
accordance  with  the  physical  nature  of  the  extra  variables.  The  initial  data  for  all  the  variables 
except  the  mass  and  volume  concentrations  and  the  old  density  are  set  to  zero.  The  initial  data 
for  these  concentrations  and  the  density  are  taken  from  the  values  determined  in  the  UINCHK 
(more  precisely,  in  the  VERCHK)  subroutine. 


4.2  Lagrangian  Block 

For  the  second  LB  part  of  the  code,  two  subroutines  of  the  main  LB  subroutine  ELSGE  are 
involved  in  the  modifications.  The  first  is  the  database  one  processing  subroutine,  ELSGD, 
and  the  second  is  the  main  processing  subroutine  ELSG  [5,  9, 10, 19], 

4.2.1  ELSG  modifications 

A  database  (DMB)  subroutine  dealing  with  the  ELSG  subroutine  is  modified  to  increase  the 
number  of  scratch  arrays  that  are  necessary  for  processing  all  the  extra  variables  by  setting  the 
relevant  variable  (the  number  of  scratch  arrays)  to  a  larger  value  for  the  present  model.  The 
ELSG  subroutine  that  is  the  main  CE  driving  part  of  the  CTH  code  processes  the  stress  tensor 
equations  taken  in  the  Jaumann  derivative  form  [20].  The  equations  for  the  deviatoric  part  of 
the  stress  tensor  are  an  analogue  of  the  equations  (2)  of  the  model.  The  initial  stage  of  the 
ELSG  modifications  extracts  both  the  symmetric  and  the  skew-symmetric  parts  of  the  velocity 
gradient  at  an  advanced  half  time  step  [20]  t  =  fn+%,  density  at  the  new  time  step  t  =  tn+ i  and 
pressure,  temperature,  and  stress  deviators  at  an  old  time  step  (at  the  start  of  the  time  cycle 
calculation)  t  =  tn.  Extra  variables,  namely,  mass,  volume,  and  entropy  disequilibrium 
concentrations,  as  well  as  'specific'  strain,  the  strain  disbalance,  and  density  and  entropy  at  the 
old  time  step  are  retrieved  and  the  extra  variable  update  is  conducted.  As  a  result,  the  main 
driving  subroutine  VERDRV  specified  for  local  DMB-extracted  arguments  as  VERSIG 
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performs  the  following:  i)  within  VERSIG,  the  extra  variables  for  mass,  volume  and  entropy 
disequilibrium  concentrations  are  calculated  by  subroutine  VEREXD,  and  ii)  the  stress 
deviator  for  the  second  constituent  is  calculated  using  equations  deduced  similarly  to  [5]  in 
the  manner  that  the  system  (21)  is  obtained  where  G2  should  be  taken  at  an  advanced  half  time 
step.  These  calculations  are  followed  by  calculations  of  the  'specific'  strain  and  strain 
disbalance  using  subroutine  VERRLX  for  calculation  of  the  stress  return  to  the  yield  surface. 

Summarizing,  the  model  modifications  of  the  ELSG  subroutine  include  the  calculation  of:  i) 
the  modulus  G2  at  t  =  fn;  ii)  the  invariant  d  of  the  tensor  af)  from  (11)  at  t  =  tn;  iii)  the  Gibbs 
potentials  and  the  potentials  A  and  Ilo  from  (8)  at  t  =  fn;  iv)  the  mass,  volume  and  entropy 
disequilibrium  concentrations  at  t  =  tn+ 1  using  the  subroutine  VEREXD;  v)  the  shear  stresses 
Sif2)  at  t  =  tn;  vi)  the  modulus  G2  at  t  =  fn+%,  as  an  average  of  Gz  at  the  old  and  new  time  steps; 
vii)  an  intermediate  stress  deviator  s’y  using  the  elasticity  equations  (ignoring  the  relaxation 
terms  in  the  stress  equations);  viii)  the  intermediate  stress  intensity  according  to  (s’)2  =  s’y  -s’y 
for  the  second  phase;  ix)  the  scaling  factor  of  the  stress  deviator  (the  yield  return)  using  the 
subroutine  VERRLX;  and  x)  the  stress  tensor  s,y(2)  at  t  =  tn+ 1  followed  by  calculations  of  e,y  and 
A  ij. 

Briefly  specifying  VEREXD,  the  mass  concentration  does  not  change  in  the  present 
implementation  due  to  absence  of  mass  exchange  between  the  phases  in  the  porous  material; 
the  volume  concentration  is  calculated  from  (4)  as  dd/dt  =  -  ip  within  the  Lagrangian  step  with 
0  taken  at  the  new  time  step  inside  ip  defined  from  (9),  (23);  and  the  entropy  disequilibrium  is 
calculated  from  (5)  as  rfj/  dt  =  -  co  within  the  Lagrangian  step  with  j  taken  at  the  new  time 
step  inside  co  defined  from  (9),  (18). 

Calculations  of  the  subroutine  VERRLX  are  practically  identical  to  the  stress  relaxation 
calculations  realised  for  the  decoupled  Maxwell- type  model  and  described  in  [5].  Namely, 
using  the  elastic  trial  stress  s’  mentioned  above  as  an  initial  data  point  within  the  Lagrangian 
step  for  the  equation  ds/dt  =  -  s/tz,  we  can  calculate  the  stress  intensity  s  at  t  =  tn+ 1  (see  [5]  for 
details).  The  stress  deviators  are  calculated  afterwards  in  the  standard  fashion: 

sn+1ij  =  ( Sn+ V  S”)  •  Snij  . 

In  the  constitutive  equation  calculations  outlined  above,  new  values  of  the  concentration 
parameters  and  the  stress  invariant  at  t  =  tn+ 1  are  calculated  by  iterations,  using  Newton's 
method. 

Thus,  the  whole  set  of  the  extra  variables  c,  0,  j,  et]  and  X„  is  calculated  in  the  present 
subroutine  for  the  Lagrangian  step.  However,  the  variables  have  to  be  updated  on  the 
Eulerian  remap  step  after  advection  and  used  afterwards  for  calculation  of  the  remaining 
thermodynamic  parameters  using  EOS  at  the  end  of  the  time  cycle  calculation. 


4.3  Eulerian  Remap  Block 

The  modification  in  the  ERB  block  deals  with  the  EREB  subroutine  [5, 9, 10, 19] .  Three  parts  of 
this  modification  deal  with  i)  preparation  of  the  extra  variables  for  the  thermodynamic 
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processing;  ii)  EOS  modification  to  the  present  model's  EOS;  and  iii)  adjustment  of  the  extra 
variables  to  the  present  thermodynamics  state  using  the  advected  values  [21]. 

4.3.1  EREB  modifications 

The  first  part  of  the  modifications  calls  subroutine  VERSWP.  The  subroutine  replaces  values 
of  the  whole  set  of  13  thermodynamic  extra  variables  defined  in  EOSVEK  and  accessed  via  the 
EOS  substitute  subroutines  by  values  of  the  first  13  (out  of  the  whole  set  of  15)  variables 
defined  in  the  UINSIV  subroutine  and  accessed  via  the  CE  subroutines.  These  variables  have 
been  updated  in  the  ELSG  subroutine  on  the  Lagrangian  step  (see  the  previous  subsection). 

The  second  part  of  the  modifications  is  dealing  with  EOS  used  in  the  following  two  forms 
admissible  by  the  CTH  code:  the  energy  form  and  the  temperature  form  (see  [22]).  The 
corresponding  substitute  subroutines  EOSVEV  and  EOSVES  have  been  modified  for  the 
present  EOS.  In  order  to  satisfy  the  variable  choice  requirement,  the  present  EOS  needs  to  be 
reformulated  in  the  density-energy  and  density-temperature  forms.  This  is  done  by  using  (6), 
(7),  (10),  (11)  and  (14).  To  facilitate  the  calculations,  entropy  is  calculated  in  these  subroutines 
from  the  internal  energy  or  temperature  by  iterations  using  Newton's  method.  The 
subroutines  must  calculate  the  following  derivatives  of  dependent  thermodynamic  variables, 
namely,  (dp/ bp)r,  (dp/ dT)p,  (be/ dp)r,  and  (be/ dT)P.  Here  energy,  e,  is  a  function  of  p  and  T.  The 
derivatives  are  obtained  by  the  standard  thermodynamics  rules  using  the  internal  energy 
potential  e(p,  S )  from  derivatives  over  p  and  S.  For  example, 

(be /dp) T  =  (be/bp)s  -  (be/ bS)P(bT/ bp)s/  (bT / bS)  P  . 

The  remaining  derivatives  are  obtained  in  a  similar  manner. 

After  the  energy  update  at  the  thermodynamics-processing  module  of  the  Eulerian  remap 
block  in  EREB,  the  last  part  of  the  modifications  finalises  preparation  of  the  extra  variables  for 
the  next  time  step.  Specifically,  using  the  internal  energy,  density  and  advected  stress  and  the 
extra  variables  that  are  common  to  the  EOS  and  EP  groups,  entropy  is  calculated  from  EOS 
and  along  with  the  density  is  stored  as  the  two  last  variables  of  the  EP  group  in  order  to  have 
these  variables  in  the  ELSG  subroutine  at  the  start  of  processing  the  next  time  cycle. 


5.  Shock  Propagation  in  Highly  Porous  Material 

The  first  set  of  test  calculations  simulates  plate  impact  experiments  in  the  one-dimensional  set¬ 
up.  This  modelling  is  aimed  to  check  out  if  the  present  implementation  is  capable  of 
describing  the  anomalous  behaviour  of  highly  porous  materials.  This  behaviour  is  well  known 
for  the  majority  of  porous  substances  at  high  and  extreme  porosities  and  manifests  itself  at  the 
shock  loading  of  porous  samples  as  decreasing  density  behind  the  shock  front  when  the 
pressure  increases  in  a  certain  pressure  range. 
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5.1  Parameters  and  Hugoniot  of  Silica 

The  specific  purpose  of  the  present  implementation  was  the  description  of  sand  subject  to 
loading  by  a  detonating  charge.  Therefore,  this  specific  porous  material  is  under  consideration 
in  the  present  work.  A  further  modification  of  the  implementation  for  the  description  of  the 
shock  behaviour  of  a  few  metallic  powders  with  the  constitutive  relations  taken  from  [4, 6-7]  is 
a  routine  procedure.  Turning  to  sand,  the  solid  constituent  of  this  porous  mixture  is  silica.  The 
pressure-density  dependences  behind  the  shock  front  (Hugoniot)  are  reported  for  many 
porous  substances  including  porous  silica.  Hugoniots  for  silica  taken  from  paper  [11]  are 
shown  in  Fig.  1  at  the  initial  material  densities  po  of  0.4  and  1.55  g/cm3.  It  is  seen  from  Fig.  1(a) 
that  at  the  high  porosity  m  more  than  6  (m  =  pos/ po)  if  po  =  0.4  g/ cm3  (for  polycrystalline  quartz 
as  the  solid  constituent,  m  =  pos/ po  =  2.65/0 .4  =  6.625),  the  Hugoniot  demonstrates  anomalous 
behaviour  in  the  pressure  range  between  5  and  20  GPa.  Besides,  the  silica  material  appearing 
in  nature  as  quartz  is  a  material  manifesting  more  complex  behaviour  than  many  metals.  In 
particular,  silica  is  subject  to  a  few  high-pressure  phase  transitions  to  polymorphs  that  have  a 
significantly  higher  reference  density  than  crystalline  or  amorphous  quartz  at  the  ambient 
pressure  (see  Fig.  1(b)  for  schematics  of  the  static  pressure-temperature  phase  diagram). 


Figure  1.  Experimental  Hugoniots  of  porous  silica  at  two  initial  densities  [11]  (a)  and  a  simplified 
schematic  of  the  silica  phase  diagram  (b) 

It  should  be  noted  that  the  simplified  diagram  in  Fig.  1(b)  ignores  numerous  high-temperature 
polymorphs  for  the  sake  of  simplicity  because  reference  densities  of  these  polymorphs  are 
quite  close  to  that  of  the  quartz  at  the  ambient  pressure  and  temperature  (the  a-quartz  for  the 
poly  crystalline  material).  Comparing  the  a-quartz  with  its  high  pressure  polymorphs,  the  first 
high  pressure  polymorph,  coesite,  shows  at  least  a  10%  reference  density  increase  and  the 
second  one,  stishovite,  more  than  60% .  Taking  into  account  a  phase  transition  within  the  solid 
constituent  would  require  a  three-phase  model  such  as  [1]  for  a  single  transition  and  a  multi¬ 
phase  model  with  four  or  more  phases  for  additional  phase  transitions.  The  present 
implementation  takes  a  simplified  approach  with  a  single  solid  phase  in  order  to  remain 
within  the  two-phase  model  framework. 

The  ability  of  the  present  model  to  describe  behaviour  of  highly  porous  metals  has  been 
demonstrated  in  papers  [4,  7],  using  an  in-house  code  based  on  the  Godunov  method.  One- 
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dimensional  calculations  with  this  code  were  verified  against  plane  impact  data  on  the  shock 
behaviour  of  a  few  porous  substances.  Because  of  the  simplified  description  of  sand  in  the 
present  work,  the  following  flowchart  will  be  used  for  verification  of  the  present 
implementation.  For  the  critical  case  of  a  highly  porous  silica,  a  one-dimensional  calculation 
will  be  conducted  using  the  code  [4, 7],  After  that,  results  of  the  calculations  will  be  compared 
with  the  modelling  results  obtained  with  CTFI  using  the  present  implementation. 

First  of  all,  the  material  constants  for  the  two-phase  model  have  to  be  specified  in  a  form  that 
is  practically  identical  for  both  the  code  [4,  7]  and  the  present  CTFI  implementation.  The 
parameters  to  be  used  are  listed  in  Subsection  4.1.3.  Amongst  them,  a  number  of  constants  are 
well  determined  such  as  the  sound  velocity  and  polytropic  exponent  of  the  air:  Co  =  0.34  km/ s 
and  y  =  1-4,  and  the  ambient  pressure,  po  =  1  bar.  For  quartz,  the  heat  capacity  is  taken  from 
[23]  cvs  =  0.75  J /  g/°K  with  the  solid  density  pos  =  2.65  g/ cm3  for  polycrystalline  quartz  and  pos 
=  2.2  g/ cm3  for  fused  quartz.  The  initial  density  of  porous  silica  to  be  used  in  the  one¬ 
dimensional  calculations  is  po  =  0.4  g/cm3.  The  shear  modulus  G  =  31.2  GPa  and  the  bulk 
modulus  K  =  37.2  GPa.  The  modulus'  pressure  derivatives  from  [24]  along  with  the  Hugoniot 
data  [25]  give:  ao  =  4.057  km/ s;  bo  =  3.43  km/ s;  ao  =  0.1;  and  /L  =  1.  While  varying  somewhat 
with  density  in  reality,  the  Griineisen  coefficient  is  taken  as  a  constant  po  =0.55  from  [26]  for 
the  pressure  range  of  interest.  The  strength  data  for  quartz  vary  in  the  literature  within  a  quite 
wide  range  from  tens  of  MPa  to  a  few  GPa.  Therefore,  as  a  plausible  choice  in  the  static- 
dynamic  range  for  a  quite  wide  span  of  strain  rates:  sx  =10-2s-1;  s2  =  103s-1,  we  take  the  yield 
limit  values  as  Yi  =  200  MPa;  Yi  =  400  MPa  at  compression.  At  tension,  the  Y-values  would  be 
more  reasonable,  in  fact,  to  relate  to  a  sandstone  rather  than  to  a  polycrystalline  or  fused 
quartz  material,  thereby  reducing  these  values  by  an  order  or  so  (e.g.  see  [27,  28]).  For  the 
shock  loading  tests  used  for  verification  of  the  behaviour  of  the  highly  porous  sand,  we  are 
mostly  interested  in  the  compression  state  of  the  material.  Therefore,  we  arbitrarily  take  Yu  = 

Yi;  Ya  =  Yi  in  the  one-dimensional  calculations.  The  solid  defects  constants  in  (15)  are  taken  as 
No  =  104  and  M  =  106.  The  compaction  kinetics  constants  in  (23)  are  chosen  to  be  identical  to 
those  from  [6]:  0oo  =  0.25;  pc  =  1.25  GPa;  nc  =  4;  ipoo  =  15  s/ m2.  For  the  heat  transfer  kinetics  we 
take  the  grain  size  as  do  =  0.1  mm  and  the  heat  conductivity  constants  for  the  air  and  quartz  as 
k\  =  0.025  W  /  m/°K  and  L  =  2W/ m/°K  with  the  dimensionless  coefficient  As  =  400  (do  and  As 
are  interrelated  when  mutually  changing;  see  [4]  for  details). 


5.2  Shock  compression  calculations  using  verified  code 

For  one-dimensional  calculations  of  the  highly  porous  silica  we  choose  the  constants  above 
with  the  initial  solid  density  corresponding  to  the  polycrystalline  quartz.  Ignoring  the  phase 
transitions  in  quartz,  we  can  evaluate  'partial'  Hugoniots  of  the  porous  material  using  the 
technique  [7]  while  keeping  various  thermodynamic  parameters  at  equilibrium. 
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Figure  2.  Calculated  PTE  Hugoniot  (solid  curve)  at  the  pressure-temperature  equilibrium  [7]  compared 
with  the  experimental  Hugoniot  data  (points  [11]) 

Traditionally,  the  Hugoniot  of  a  multi-phase  material  is  a  locus  of  states  behind  the  shock 
front  that  are  considered  to  be  in  thermodynamic  equilibrium  between  the  material  phases. 
However,  as  noted  in  [7],  this  is  not  always  the  case.  At  certain  levels  of  loads  some 
thermodynamic  parameters  between  the  phases  may  be  in  equilibrium,  whereas  the  others  are 
not.  Thus,  a  variety  of  Hugoniots  can  be  derived  and  depending  on  the  meso-mechanics  of  a 
multi-phase  material  one  Hugoniot  from  the  variety  can  supersede  another,  resulting  in  a 
composite  Hugoniot  for  a  real  material  under  shock  loading  (see  [1,  4,  6,  7]  for  details  and 
examples).  Specifically,  for  some  typical  thermodynamic  parameters  Hugoniots  are  referred  to 
as  PE  (Pressure  Equilibrium)  and  PTE  (Pressure-Temperature  Equilibrium)  Hugoniots  (other 
Hugoniots  are  also  possible,  see  [1, 4, 7]  for  details).  For  the  present  case,  where  the  nature  of 
equilibrium  is  changing  during  shock  loading,  the  anomalous  pressure-density  behaviour 
behind  the  shock  front  (the  PTE  Hugoniots  are  usually  anomalous  for  highly  porous  mixtures) 
supersedes  the  conventional  behaviour  (the  compaction  stage).  For  illustration,  the  calculated 
analytical  PTE  Hugoniot  of  sand  at  po  =  0.4  g/ cm3  shown  in  Fig.  2  by  a  solid  curve  manifests 
the  anomalous  behaviour.  However,  as  mentioned  above,  this  behaviour  is  only  a  part  of  the 
whole  picture.  Different  porous  mixtures  realize  different  mechanisms  of  transition  from  one 
regime  to  another. 


The  highly  porous  silica  demonstrates  a  relatively  low  point  of  transition  to  the  anomalous 
path.  This  is  likely  to  happen  immediately  after  compaction.  Thus,  the  experimental  points  in 
Fig.  2  demonstrate  the  anomalous  behaviour  above  5  GPa  and  conventional  behaviour  below 
this  point,  where  the  compaction  curve  reflects  the  conventional  relationship  of  density  rise 
with  pressure  increase.  It  should  be  kept  in  mind,  that  the  density  shift  of  the  calculated 
Hugoniot  from  the  experimental  points  in  Fig.  2  is  possibly  associated  with  a  phase  transition, 
and  it  is  not  taken  into  account  for  the  sake  of  simplicity.  However,  an  attempt  to  simulate  the 
switch  from  the  conventional  to  the  anomalous  behaviour  observed  in  the  experiments  shown 
in  Fig.  2  will  be  undertaken. 
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Figure  3.  Numerical  results  for  the  porous  sand  using  the  code  [4] 

The  first  calculation  is  conducted  with  the  one-dimensional  in-house  code  [4,  7]  based  on  the 
Godunov  method  [29] .  This  code  has  been  verified  for  a  few  porous  and  condensed  two-phase 
mixtures  [4,  6,  7]  using  relations  including,  but  not  limited  to  the  constitutive  functions 
described  above.  Numerical  results  obtained  with  this  code  are  shown  in  Fig.  3  for  the  high 
velocity  impact  set-up  with  a  copper  flyer  plate  and  a  target  representing  a  layer  of  sand. 
Thicknesses  of  the  flyer  plate  (hj )  and  the  sand  sample  (/is)  are  h/=  2.5  cm,  hs  =  3.5  cm  ( a-b )  and 
h/=  1.5  cm,  hs  =  2.5  cm  (c-f).  The  velocities  of  impact,  Uo,  in  km/ s  for  Figs.  3(a-f)  are  0.5  (a);  1 
(b);  3  (c);  5  (d);  7  (e),  and  9  (/).  The  calculation  is  conducted  in  adaptive  grids  with  the  combined 
computational  domain  (flyer  plate  plus  target)  rescaled  back  for  the  sake  of  visualisation  to 
the  initial  domain  as  shown  in  Fig.  3.  The  free  and  contact  interfaces  in  these  calculations  are 
Lagrangian  and  internal  nodes  -  Eulerian  (see  [29]  for  details).  For  simulation  of  the  flyer 
plate,  the  viscoelastic  model  [17]  has  been  used  with  the  material  constants  for  copper  taken 
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from  [13].  Snapshots  of  the  density  and  pressure  profiles  are  shown  at  the  same  time  points 
for  each  of  the  impact  velocities  (the  time  marks  are  shown  above  each  of  the  graphs).  Tracing 
the  pressure  and  density  states  behind  the  shock  front  in  sand,  it  is  seen  that  when  pressure  is 
rising,  density  is  first  increasing  in  the  pressure  range  below  5  GPa  (plots  in  Fig.,  3 (a-c))  and 
then  decreasing  at  further  increase  of  pressure  (plots  in  Fig.,  3 (d-f)).  This  behaviour  agrees 
quite  well  with  the  experimental  behaviour  shown  in  Fig.  2,  following  the  trend  of  the  initial 
conventional  compaction  path  superseded  by  the  anomalous  PTE  Flugoniot. 


5.3  Shock  compression  calculations  using  the  CTH  implementation 

Calculations  with  CTH  using  the  present  implementation  require  specification  of  the  sand 
material  in  the  both  sections  of  the  input  deck  associated  with  the  EP  and  EOS  blocks.  The 
corresponding  specification  in  the  'eos'  section  refers  to  VE  EOS  subroutines  that  are 
substituted  in  the  present  implementation  by  the  subroutines  described  in  the  previous 
section  of  the  report. 
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Figure  4.  Comparison  of  the  CTH  calculations  in  the  underlying  (a)  and  control  (b)  set-ups  with  the 
profiles  calculated  using  the  code  [4,  7]  (c) 


This  reference  forces  CTH  to  retrieve  information  about  any  material  from  the  EOS_data 
material  database  associated  with  the  VE  model  (for  example,  in  a  fragment  of  the  input  deck 
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below,  the  'pmma'  material  is  referred  to).  Because  of  the  substitution  this  information  cannot 
actually  be  utilized  when  the  implementation  is  active  with  the  exception  of  the  first  three 
constants.  In  fact,  only  these  three  constants  are  important  for  the  corresponding 
implementation  subroutines  in  order  to  start  proper  processing  of  the  EOS  input  and 
initialisation  blocks.  After  processing  of  the  EP  input  block,  the  whole  set  of  constants  is 
determined  and  correctly  placed  for  subsequent  calls  of  the  EOS  modules.  These  three 
constants  taking  the  roles  of  the  mixture  density,  and  densities  of  the  gaseous  and  solid  phases 
are  'RO',  'TO',  and  'CV',  respectively,  referring  to  the  VE  input  nomenclature.  The  processing 
of  the  EP  input  data  is  conducted  in  quite  a  conventional  way  in  CTH  through  reference  to  the 
present  model  name,  'VER'.  Specification  of  one  of  the  constants,  namely,  the  material  density 
(unless  this  density  is  coincident  with  the  one  from  the  dataset  called  by  the  material  name 
specification),  is  mandatory,  whereas  the  whole  set  of  the  EP  input  constants  is  summarised  in 
the  VP_data  file.  This  data  set  is  taken  from  the  rows  of  the  VP_data  file  identified  by  the 
material  name  and  it  contains  37  constants  described  in  Subsection  4.1.3. 

For  example,  assuming  that  the  material's  '2'  name  is  'SAND',  and  its  initial  density  po  =  0.4 
g/ cm3,  the  required  fragment  of  material  specification  in  the  input  deck  for  the  material 
described  by  the  present  two-phase  model  should  appear  as  follows 

eos 

mat2  ve  pmma  R0=0.4  T0=0 . 00121107266  CV=2.65 
endeos 

epdata 

matep  2  VER  SAND  RHO  =  0.4 
endepdata 

To  remind,  the  initial  density  of  the  air  pog  =  0.00121107266  g/ cm3  is  actually  calculated  on  the 
processing  stage  of  the  EP  input  block  and  specified  here  only  for  the  correct  initial  processing 
of  the  EOS  input  and  initialisation  blocks.  Therefore,  this  density  is  calculated  beforehand  and 
used  as  the  datum  'TO'  above.  For  the  same  purpose,  the  initial  density  of  the  solid  phase  pos  = 
2.65  g/ cm3  is  also  taken  from  the  datum  'CV'  above.  It  should  be  noted  that  this  twist  is 
unnecessary  if  proper  access  to  the  EOS  block  is  available.  In  that  case,  new  EOS  name  can  be 
specified  and  used  for  the  present  model.  Therefore,  the  model  could  be  realised  with  its  own 
specifications  instead  of  using  the  name,  VE,  of  the  subroutines  substituted  and  the  names  of 
VE  variables  for  substitution  in  the  data  set,  when  the  relevant  'eos'  section  in  the  input  deck 
is  activated. 

The  high-velocity  plane  impact  set-up  for  CTTI  employs  similar  input  conditions  as  in  the 
previous  calculation.  The  computational  one-dimensional  Eulerian  domain  is  separated  into 
two  zones:  i)  a  flyer  plate  with  the  thickness  hf  of  either  2.5  or  3.5  cm  depending  on  the  impact 
velocity,  similarly  to  the  previous  subsection  set-ups;  and  ii)  the  sand  sample  zone  (a  target) 
with  the  thickness  hs  of  either  1.5  or  2.5  cm,  respectively.  Due  to  the  purely  Eulerian  nature  of 
the  code,  the  boundary  conditions  require  more  attention  than  in  the  previous  case.  In  the 
CTH  simulations  below,  the  left  boundary  condition  was  usually  taken  as  the  inflow  condition 
that  enabled  the  rarefaction  wave  from  the  rear  side  of  the  flyer  plate  to  come  to  the  projectile- 
target  interface  later  than  in  the  case  of  the  free  surface  condition.  However,  this  choice  did 
not  affect  the  density  value  behind  the  shock  front.  This  point  has  been  checked  out  by  special 
control  calculations  with  the  free  surface  conditions. 
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Figure  5.  CTH  calculations  of  the  density  profiles  in  the  porous  sand  at  various  impact  velocities 

In  the  CTH  calculations,  a  model  for  the  description  of  the  copper  material  of  a  flyer  plate  is 
the  Steinberg-Guinan-Lund  model  [3].  Results  of  comparison  of  the  CTH  calculations  using 
the  present  implementation  at  Uo  =  1  km/  s  with  the  set-up  described  above  (a),  with  a  control 
set-up  with  the  free  surface  boundary  conditions  (b),  and  the  calculation  from  the  previous 
subsection  (c)  are  shown  in  Fig.  4.  It  is  seen  that  all  these  calculations  result  in  nearly  identical 
density  profiles  in  the  sand  sample.  The  pressure  profiles  for  the  same  boundary  conditions 
(Figs.  4 (b)  and  (c))  are  close.  It  should  be  kept  in  mind  the  visualisation  convention  in  Fig.  4(c), 
where  the  external  free  boundaries  are  rescaled  back  to  the  initial  positions  with  the  contact 
interface  having  changed  its  position  between  the  external  free  boundaries  by  t  =  12  psec, 
whereas  the  CTH  visualisation  in  Figs.  4 (a-b)  refers  to  the  actual  boundary  positions  within 
the  fixed  grid.  The  observed  pressure  jump  at  the  projectile- target  interface  is  associated  with 
the  strength  effects  because  the  stress  equilibrium  at  the  interface  is  satisfied  only  for  the 
longitudinal  component  of  stress. 

The  density  profiles  from  the  CTH  calculations  using  the  present  model  implementation  at  the 
same  impact  velocities  as  in  the  previous  subsection  (Uo  in  km/ s  are  0.5(a);  1(b) ;  3(c);  5(d);  7(e); 
and  9(f))  are  summarized  in  Fig.  5  for  the  boundary  conditions  of  the  underlying  set-up  (as 
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that  of  Fig.  4(a)).  When  comparing  the  corresponding  profiles  from  Figs.  3  and  5,  it  is  seen  that 
the  profile  features  in  sand  such  as  the  peak  density  and  the  character  of  ramping  are 
practically  identical  for  these  two  series  of  calculations.  Thus,  this  comparison  gives  us  a 
confidence  in  the  adequacy  of  the  present  numerical  realisation  within  the  CTFI 
implementation  framework. 


Figure  6.  Summary  of  the  Hugoniot  densities  calculated  with  CTH  using  the  present  model 
implementation 


The  numerical  data  summarised  from  Fig.  5  are  drawn  in  Fig.  6  as  a  single  graph  with  the 
straight  lines  characterising  the  trend  of  the  behind-shock  density  change  with  the  pressure 
rise  from  left  to  right.  This  trend  shows  the  physical  adequacy  of  the  implementation  (within 
the  accuracy  of  chosen  EOS),  when  comparing  the  calculated  behind-shock  (Hugoniot) 
densities,  following  a  composite  Hugoniot  from  the  compaction  pressure-density  curve 
superseded  by  the  anomalous  PTE  Hugoniot,  which  correlates  quite  well  with  the  Hugoniot 
response  observed  experimentally  [11]  (see  Fig.  2). 


5.4  Shock  compression  calculations  using  models  available  in  CTH 

The  CTH  calculations  using  the  present  model  can  also  be  compared  with  CTH  calculations 
employing  porous  material  models  available  in  the  CTH  model  database  [3].  The  CTH  models 
that  may  associate  most  closely  the  thermodynamic  states  of  a  material  with  EOS 
experimental  data  are  those  of  the  SESAME  family  [30].  The  SESAME  EOS  database  describes 
the  material  state  response  in  a  tabulated  form  with  necessary  interpolations  between  the  state 
database  nodes.  However,  a  disadvantage  of  this  approach  is  that  any  out-of-node  state  is  a 
calculated  extrapolation  of  the  tabulated  data  and,  if  not  in  close  proximity  to  the  table  data, 
this  approximation  relies  on  extrapolation  algorithms  rather  than  on  physics  laws.  A  relevant 
model  from  the  SESAME  database  is  the  DRY  SAND  model  (in  fact,  a  table  extrapolation) 
complemented  with  the  P-alpha  EOS  reduction  rule  [3]  used  for  the  porous  material 
description.  Because  this  description  extrapolates  a  behaviour  between  the  EOS  data  nodes 
tabulated  by  experiment,  it  is  expected  that  the  pressure  density  abnormalities  are  somewhat 
covered  by  the  DRY  SAND  model  prediction.  The  second  model  is  the  P-X  model  by  Grady 
that  describes  the  non-linear  response  using  a  composite  EOS  in  a  kinetic  manner  managed  by 
the  compaction  parameter  X  [31]. 
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Figure  7.  Comparison  of  the  numerical  results  for  the  present  model  (a)  with  those  for  the  DRY  SAND 
(b)  and  P- X  (c)  models  ofCTH 


The  pressure  and  density  profiles  taken  from  the  previous  subsection  and  calculated  using  the 
present  two-phase  model  at  the  impact  velocity  of  Uo  =  1  km/ s  are  shown  in  Fig.  7(a).  These 
calculations  are  compared  with  the  profiles  within  the  same  set-up,  which  are  calculated  using 
the  DRY  SAND  model  (Fig.  7(b))  and  P-X  model  (Fig.  7(c)).  It  is  seen  that  the  tabular  EOS 
response  based  on  the  SESAME  approach  in  Fig.  7(b)  manifests  an  oscillating  pressure  profile 
behind  the  shock  front  fed  into  the  pressure  pulse  propagating  back  to  the  flyer  plate.  In  turn, 
the  response  obtained  with  the  physics-based  P-X  model  in  Fig.  7(c)  results  in  a  quite  similar 
pressure  profile  in  the  flyer  plate  when  compared  with  that  of  the  present  two-phase  model  in 
Fig.  7(a).  The  pressure  amplitude,  however,  is  obtained  from  the  model-specific  EOSs, 
therefore,  it  is  not  expected  to  be  identical  for  these  models.  At  the  same  time,  the  density 
profiles  in  sand  are  quite  different  for  all  the  three  models. 
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Figure  8.  Trend  of  the  Hugoniot  densities  versus  pressure  calculated  with  the  DRY  SAND  model 

The  density  profiles  calculated  with  CTH  using  the  DRY  SAND  model  at  the  same  impact 
velocities  as  in  the  previous  subsections  are  summarised  in  Fig.  8  in  the  same  manner  as  in 
Fig.  6.  It  can  be  seen  that  the  tabular  character  of  the  SESAME  EOS,  resulting  in  a  somewhat 
oscillating  response,  transmits  this  response  to  the  density  trend  at  high  pressure.  However,  in 
the  middle  of  the  pressure  range  the  correct  abnormality  trend  still  remains,  namely,  the 
maximum  density  is  reached  at  the  impact  velocity  of  Uo  =  3  km/ s  (the  third  graph  from  the 
left  in  Fig.  8),  corresponding  to  the  peak  pressure  of  approximately  5  GPa  (see  Hugoniot  data 
in  Fig.  2). 


X  (cm)  X  (cm)  X(cm)  X(cm)  X(cm)  X(cm) 


Figure  9.  Trend  of  the  Hugoniot  densities  versus  pressure  calculated  with  the  P- X  model 

The  next  model  prediction  with  a  regular  EOS  such  as  that  of  the  P-A  model  does  not 
demonstrate  any  abnormality.  The  corresponding  CTH  results  are  shown  in  Fig.  9  and  the 
density  trend  is  quite  conventional  with  the  density  and  pressure  increasing  concordantly. 
Obviously,  this  trend  does  not  correlate  with  the  experimental  data  shown  in  Fig.  2. 


6.  Buried  Charge  Modelling 

In  the  mine  blast  experiments,  total  momentum  is  often  recorded  by  using  a  momentum¬ 
measuring  device  (e.g.  pendulum).  However,  the  major  contribution  to  the  total  momentum  is 
due  to  loading  by  the  detonation  products,  while  the  major  contribution  to  the  stress  pulse 
transferred  to  the  target  is  due  to  the  ejecta  impact.  Therefore,  in  order  to  take  into  account  the 
highly  dangerous,  although  short-time  ejecta  contribution,  the  soil  ejecta  formation  and 
momentum  transfer  by  the  ejecta  to  a  target  should  be  considered.  As  mentioned  in  the 
Introduction,  the  velocity-density  distribution  through  the  ejecta  thickness  is  the  decisive 
information  for  evaluation  of  this  contribution.  This  information,  however,  is  not  readily 
available.  The  peak  stress  in  the  stress  pulse  travelling  through  the  target  is  associated  with 
the  impact  velocity  of  the  soil  ejecta.  The  length  of  the  pulse  and  the  pulse  fluctuations  could 
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be  obtained  from  the  velocity-density  information.  Nevertheless,  as  a  first  step  in  the 
evaluation,  the  velocity  can  accurately  be  estimated  from  high-speed  video  or  pulse  x-ray 
shots.  CTH  calculations  in  the  present  section  are  aimed  at  analysis  of  the  experimental 
observations  made  with  the  pulse  x-ray  experiments  [12]  on  the  explosion  of  lOOg  C4  charges 
buried  under  a  layer  of  sand. 


6.1  The  problem  statement 

The  experimental  results  [12]  are  simulated  within  the  axi-symmetrical  two-dimensional  set¬ 
up.  The  experiments  [12]  deal  with  a  cylindrical  volume  of  sand  of  approximately  88  cm  in 
diameter  and  70  cm  in  height  that  is  not  confined  at  the  upper  side  of  the  volume  (i.e.  it  is  in 
the  contact  with  air)  and  is  contained  laterally  and  at  the  bottom  by  a  steel  case.  The  present 
numerical  set-up  considers  one  half  of  the  whole  experimental  set-up  in  the  radial  direction 
away  from  the  symmetry  axis  with  rigid  boundary  conditions  sidewise  and  at  the  bottom. 
20  cm  of  air  space  is  within  the  chosen  computational  domain.  The  observation  area  subject  to 
graphical  output  is  40  x  40  cm  with  the  centre  at  the  intersection  of  the  symmetry  axis  and  the 
air-sand  free  surface  interface.  This  observation  area  is  marked  by  the  square  in  the  set-up 
graph  of  Fig.  10  for  the  whole  set-up  geometry  (the  actual  observation  area  in  the  calculation 
set-up  is  20  x  40  cm),  which  approximately  corresponds  to  the  experimental  observation  area 
[12]. 

Two  Depths  Of  Burial  (DOB)  for  the  High  Explosive  (HE)  charges  are  considered  in 
agreement  with  the  buried  charge  set-ups  of  [12]:  DOB  =  3  cm  and  DOB  =  8  cm.  The  C4  charge 
was  simulated  using  the  HE  burn  model  [3]  under  JWL  EOS  for  Composition  C-4  (see  [3]  for 
details)  with  the  JWL  database  choice  of  constants.  The  charge  was  selected  to  be  cylindrical 
with  a  diameter  of  6.4  cm  and  a  height  of  2  cm  [12].  Two  charge  initiation  configurations  have 
been  considered  in  [12]  with  detonation  either  (i)  at  the  top  or  (ii)  at  the  bottom  of  the  charge. 
As  found  in  [12],  the  first  option  generates  a  venting  effect  at  the  soil  interface  rendering  this 
set-up  inconclusive,  whereas  the  latter  is  free  of  such  limitations.  Therefore,  only  the  charge 
initiation  configuration  (ii)  was  considered  in  the  present  numerical  set-ups.  To  approximate 
the  experimental  set-up  the  charge  was  encased  in  a  1  mm  thick  Nylon  case  described  with  a 
standard  elasto-plastic  model  (von  Mises  criterion  is  used  for  the  plastic  flow  simulation)  with 
the  yield  limit  of  500  MPa  and  the  'NYLON'  SESAME  EOS  [3].  Similarly  to  the  experimental 
set-ups,  the  charge  was  boosted  by  a  5  x  10  mm  HE  cylinder  detonated  from  its  lowest  point  at 
the  axis  of  symmetry  in  the  numerical  set-up. 

Material  constants  for  the  sand  material  modelled  by  the  present  CTH  implementation  are  the 
same  as  those  from  the  preceding  subsection  with  the  only  adjustment  to  the  fused  quartz 
density  value  of  pos  =  2.2  g/ cm3.  This  type  of  solid  phase  for  calculations  in  this  section  is 
chosen  because  the  EOS  in  the  present  case  should  cover  a  wider  range  of  pressures,  including 
those  occurring  at  the  rarefaction.  Therefore,  the  Hugoniot  shift  to  the  higher  density  area 
seen  in  Fig.  2  might  not  be  as  critical  for  this  porous  material  with  a  relatively  low  porosity 
opposed  to  the  high  porosity  sand  analysed  at  compression  in  the  previous  section.  The  choice 
of  yield  limits  Yi  and  Yo  is  somewhat  artificial  because  the  yield  limits  are  related  only  to  the 
solid  constituent  and  are  not  generally  associated  with  the  porous  mixture.  This  choice  was 
quite  obvious  in  the  preceding  section  because  only  the  shock  compression  behaviour  was 
considered  and  the  yield  limits  were  associated  with  the  deformation  (crush)  of  solid  grains. 
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The  material  densification  can  potentially  be  described  by  the  compaction  kinetic  fitted  at 
different  modes  of  loading.  Namely,  in  the  case  of  cyclic  loading  and  unloading  of  the  porous 
material,  the  compaction  kinetic,  ip,  should  be  fitted  by  an  analysis  of  solutions  of  (21)  for 
specifically  designed  loading-unloading  tests.  Thus,  this  kinetic  plays  a  key  role  in  the 
analysis  of  this  sort  of  problem.  In  the  present  work,  we  are  forced  to  use  the  compaction 
kinetic  fitted  only  to  the  compression  mode  of  loading. 


Mesh  at  0.00e+00  seconds 


X  (cm) 

Figure  10.  Schematic  of  the  Buried  Charge  set-up 

Therefore,  we  are  trying  to  compensate  the  gap  in  determination  of  the  compaction  kinetic  for 
this  material,  which  significantly  contributes  to  the  behaviour  of  the  soil  ejecta,  by  a  variation 
of  the  yield  limits.  In  doing  so,  we  assume  that  the  behaviour  of  the  porous  material  at 
unloading  can  be  compensated  by  the  deformation  of  the  solid  constituent,  which  is  a  very 
rough  assumption.  Several  CTF1  runs  have  been  conducted  with  various  yield  limit  values  at 
compression  and  tension,  which  are  listed  in  Table  1.  For  completeness,  even  though  they  are 
unrealistic,  the  set-ups  with  a  high  strength  at  tension  are  included  in  the  Table  as  well. 

The  initial  average  density  of  sand  for  all  the  cases  calculated  in  this  section  is  pa  =  1.57  g /  cm3, 
that  corresponds  to  an  average  initial  soil  density  in  the  tests  [12],  The  same  porous  material 
models  from  the  CTF1  material  model  database  as  in  the  previous  section  have  been  selected 
for  comparative  calculations  and  listed  in  Table  1.  As  a  reminder,  the  models  are:  the  DRY 
SAND  model  from  the  SESAME  EOS  family  with  P-alpha  model  for  the  porosity  simulation 
and  the  P-A  model  by  D.  Grady  [31]. 
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Table  1.  Modelling  set-ups  and  yield  limit  values  for  the  implemented  model 


N 

Compression 

Tension  (p  <  0) 

Yi,  MPa 
(£,  =10-2s-!) 

Y2/  MPa 
(^2=103s-!) 

Yit/  MPa 
(s,  M0-V1) 

Y2t,  MPa 
(^2=103s-!) 

DOB  =  8  cm 

Runl 

37.5 

150 

37.5 

150 

Run2 

10 

40 

10 

40 

Run3 

37.5 

150 

10 

40 

Run4 

DRY  SA] 

ND  model  [3] 

Run5 

P-Lambda  model  [3] 

DOB  =  3  cm 

Run6 

37.5 

150 

37.5 

150 

Run7 

10 

40 

10 

40 

Run8 

37.5 

150 

10 

40 

Run9 

DRY  SAl 

ND  model  [3] 

RunlO 

P-Lambda  model  [3] 

Von  Mises  criterion  of  the  traditional  elasto-plastic  model  has  been  used  for  the  DRY  SAND 
runs  with  the  yield  limit  of  10  MPa  and  for  the  P-A  model  with  the  yield  limit  of  500  MPa. 
These  yield  limit  values  provide  the  best  fit  between  calculation  results  and  the  present 
experimental  data  within  a  reasonable  yield  limit  range  for  quartz. 


Table  2.  Summary  of  the  experimental  data  [12] 


N 

Time 

(ps) 

DOB 

(cm) 

Base  width 
(cm) 

1/3  width 
(cm) 

2/3  width 
(cm) 

Bubble  height 
(cm) 

12-1 

230.9 

8 

19.9 

13.0 

9.4 

2.9 

12-2 

301.1 

8 

22.9 

14.5 

9.8 

4.3 

13-1 

351.0 

8 

27.9 

16.1 

11.1 

5.2 

13-2 

401.1 

8 

29.1 

16.7 

11.9 

6.2 

14-1 

451.0 

8 

32.8 

17.7 

11.8 

6.9 

14-2 

526.1 

8 

31.1 

17.6 

12.4 

8.2 

16-1 

100.8 

3 

17.3 

11.5 

7.8 

6.5 

15-1 

100.9 

3 

17.7 

12.3 

9.8 

8.4 

17-1 

125.9 

3 

19.9 

13.1 

10.1 

8.8 

16-2 

201.1 

3 

22.0 

17.3 

13.8 

16.4 

17-2 

201.1 

3 

24.6 

17.3 

14.3 

16.5 

15-2 

201.9 

3 

24.2 

20.7 

18.5 

19.2 

The  experimental  results  [12]  are  summarised  in  Table  2  showing  the  basic  shape 
characteristics  of  sand  ejecta  for  the  two  series  of  tests  with  DOB  =  8  cm  (3  tests  with  2 
snapshots  for  each  test)  and  DOB  =  3  cm  (also  3  tests  with  2  snapshots  per  test).  The 
geometrical  parameters  taken  from  the  x-ray  shots  are  somewhat  subjective  but  the  Table  is 
useful  to  see  the  test  consistency  and  to  evaluate  the  general  appearance  and  velocity  of  the 
soil  ejecta. 
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6.2  Calculations  for  Large  Depth  of  Burial 

The  first  series  of  calculations  (runs  from  1  till  5  in  Table  1  that  are  further  referred  to  as  Runl- 
Run5)  deals  with  identical  mesh  set-ups  for  the  whole  group  of  the  five  with  DOB  =  8  cm.  This 
depth  of  burial  is  characteristic  for  conditions  of  the  camouflet  blast  where  the  blast  products 
are  contained  underground  during  the  ejecta  formation. 


Mesh  2.31c-04  seconds  Mesh  at  3.01  c-04  seconds  Mesh  at  3.5 le-04  seconds 


-20  -10  0  10  20  -20  -10  0  10  20  -20  -10  0  10  20 
X(cm)  X  (cm)  X(cm) 

Figure  11.  CTH  calculation  of  the  sand  ejecta  using  the  present  model  with  parameters  for  Runl 

Numerical  results  of  the  first  calculation  from  this  series  (see  Table  1),  Runl,  are  shown  in  Fig. 
11  at  times  approximately  corresponding  to  those  from  the  experimental  observations  shown 
in  Table  2.  It  is  seen  that  the  material  strength  of  the  solid  constituent  equally  high  at  tension 
and  compression,  which  is  approximately  correspondent  to  that  of  non-porous  quartz,  has  a 
significant  and  unrealistic  constraining  effect  on  the  soil  expansion.  It  is  clear  that  the  strength 
properties  at  tension  should,  at  least,  somehow  differ  from  those  at  compression. 

As  mentioned  in  the  earlier  comments  on  Table  1,  the  compaction  kinetic  fit  should,  in  fact,  be 
responsible  for  the  difference  in  the  material  compressibility  of  the  grained  porous  sand  at 
compression  and  at  tension.  Because  the  compaction  kinetic  selected  for  the  compression  case 
was  not  adjusted  to  the  tension  conditions  in  this  work,  the  next  step  in  the  material 
properties  adjustment  was  reduction  of  the  material  strength  for  the  solid  constituent. 
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Mesh  a i  2.31e-04  seconds  Mesh  at  3  Ole-CM  seconds  Mesh  at  3  51e-04  seconds 


-20  -10  0  10  20  -20  -10  0  10  20  -20  -10  0  10  20 
X(cm)  X  (cm)  X(cm) 

Figure  12.  CTH  calculation  of  the  sand  ejecta  using  the  present  model  with  parameters  for  Run2 

For  analysis  of  the  influence  of  the  strength  reduction,  the  next  calculation  from  Table  1,  Run2, 
is  conducted  with  the  conditions  of  low  yield  limit  at  both  compression  and  tension.  From  the 
results  of  this  run,  we  can  observe  an  extension  of  the  soil  ejecta  both  at  the  ejecta  base  and 
upwards  resulting  in  a  significant  increase  of  the  bubble  height  when  compared  with  the 
previous  set-up.  Results  of  this  modelling  using  the  present  implementation  are  summarised 
in  Fig.  12,  similarly  to  Runl  in  Fig.  11. 

For  completeness  of  choice  of  the  yield  limit  variations,  the  last  option  calculated  within  this 
series  is  calculation  Run3  (Table  1)  with  high  yield  limits  for  the  solid  constituent  at 
compression  corresponding  to  the  set-up  of  Runl  and  the  low  yield  limits  at  tension 
corresponding  to  the  set-up  of  Run2.  Numerical  results  for  this  run  are  summarised  in  Fig.  13 
similarly  to  the  previous  cases. 
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Mesh  at  2  3 lc-04  seconds  Mesh  at  3  0 Ic-04  seconds  Mesh  af  3.5 le-04  seconds 


-20  -10  0  10  20  -20  -10  0  10  20  -20  -10  0  10  20 
X(cn»)  X(cm)  X(cm) 

Figure  13.  CTH  calculation  of  the  sand  ejecta  using  the  present  model  with  parameters  for  Run3 

Wrap  up  calculations  of  this  series,  Run4  and  Run5,  have  been  conducted  with  CTH  using  the 
material  database  models  employing  the  DRY  SAND  model  based  on  the  SESAME  EOS  and 
the  P-A  model.  Corresponding  numerical  results  are  summarised  in  Figs.  14  and  15, 
respectively.  The  P-A  model  results  describe  shapes  of  the  ejecta  bubble  quite  poorly  and, 
therefore,  only  3  snapshots  are  shown  in  Fig.  15  for  comparison  with  the  numerical  results 
obtained  with  the  other  models. 

Analysing  the  CTH  results  obtained  with  the  present  implementation  (Figs.  11-13),  we  can  see 
that  the  case  of  Runl  with  very  high  constraining  effects  in  Fig.  11  and  the  case  of  Run2  with 
significant  relieving  effects  in  Fig.  12  are  separated  by  the  case  for  Run3.  When  comparing  the 
Run3  case  (Fig.  13)  with  that  of  Runl  (Fig.  11),  it  is  seen  that  the  high  yield  limits  at  tension 
(with  the  same  yield  limits  at  compression)  significantly  affect  the  bubble  height.  At  the  same 
time,  when  comparing  the  Run3  case  (Fig.  13)  with  Run2  (Fig.  12),  it  is  seen  that  with  the 
increase  of  the  yield  limit  to  the  high  values  at  compression  (with  the  same  yield  limits  at 
tension),  the  bubble  height  is  almost  unaffected  and  the  base  size  is  slightly  restrained. 
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Mesh  al  2.31c-04  seconds  Mesh  at  3.01c-04  seconds  Mesh  al  4  0  le-04  seconds 


X  (cm)  X  (an)  X  (cm) 

Figure  14.  CTH  calculation  of  the  sand  ejecta  using  the  DRY  SAND  model  (Run4) 


Figure  15.  CTH  calculation  of  the  sand  ejecta  using  the  P- X  model  (Run5) 


In  order  to  compare  the  results  of  the  present  series  with  the  experiments,  we  superimpose  the 
calculated  external  ejecta  contours  from  Figs.  11-15  with  the  experimental  snapshots  [12]  at  the 
late  stages  of  the  ejecta  formation.  The  results  for  the  whole  set  of  runs,  Runl-Run5  (Table  1), 
are  summarised  in  Fig.  16  (corresponding  contours  of  the  CTH  calculations  are  referred  to  by 
the  run  names  in  each  frame  of  Fig.  16).  From  the  comparison,  it  is  seen  that  the  numerical 
results  for  Run2,  Run3,  and  Run4  are  in  reasonable  agreement  with  the  experiment.  The 
constraining  effect  seen  in  the  Runl  graphs  seems  to  be  unrealistic,  and  choice  between  Run2 
and  Run3  can  be  done  only  after  the  appropriate  fitting  of  the  compaction  kinetic  which 
would  also  describe  the  unloading  behaviour  of  the  material. 
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Shot  14-1.  451  ns  8  cm  DoB 

Figure  16.  Comparison  of  experimental  snapshots  at  late  stage  with  the  CTH  results 

The  results  of  Run2  are  drawn  versus  all  available  experimental  snapshots  in  Fig.  17.  It  can  be 
seen  that  the  CTH  calculation,  employing  the  present  implementation  and  selected 
parameters  for  the  Run2  set-up,  provides  very  good  agreement  with  the  experiment.  It  is 
interesting  to  note  that  all  available  model  database  calculations  provide  higher  ejecta  velocity 
than  that  observed  in  the  camouflet  blast  experiments  (with  a  fairly  small  excess,  though,  for 
the  DRY  SAND  model).  At  the  same  time,  CTH  calculations  (not  described  here)  using  the 
present  implementation  with  even  lower  yield  limits  do  not  result  in  a  further  increase  of  the 
ejecta  velocity. 


Shot  14-2.  526. 1  )is  8  cm  DoB 


Shot  14-2.  526. 1  ns  8  cm  DoB 


Shot  14-2.  526. 1  (is  8  cm  DoB 


526. 1  (is  8  cm  DoB 


Run5 


Shot  14-2. 
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Shot  13-1.  351  (is  S  cm  DoB  Shot  13-2.  401.1  (is  S  cm  DoB 


Shot  14-1.  45 1  |is  8  cm  DoB  Shot  14-2.  526. 1  jis  8  cm  DoB 


Figure  17.  Comparison  of  the  Run 2  calculations  with  the  experimental  snapshots  [12] 


6.3  Calculations  for  Small  Depth  of  Burial 

When  depth  of  burial  decreases,  the  soil  expansion  regime  changes  from  a  confined 
underground  blast  (the  camouflet  blast)  to  soil  ejection  accompanied  by  an  outburst  of  the 
detonation  products  through  the  soil  surface.  However,  when  the  impact  effect  of  the  ejecta  on 
a  target  is  significant,  the  particles  and  detonation  products  are  still  relatively  close  to  each 
other.  Therefore,  velocity  equilibrium  between  products  and  particles  can  be  assumed  for  this 
case,  which  allows  us  to  describe  the  ejecta  behaviour  using  the  same  models  as  in  the 
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previous  subsection.  The  outburst  regime  realised  at  DOB  =  3  cm  corresponds  to  the  next 
characteristic  mesh  set-up  selected  for  the  second  series  of  CTH  calculations. 


Mesh  at  1  01c-04  seconds  Mesh  at  I  26c-04  seconds  Mesh  al  2  02e-04  seconds 


-20  -10  0  10  20  -20  -10  0  10  20  -20  -10  0  10  20 
X(cm)  X  (cm)  X(cm) 

Figure  18.  CTH  calculation  of  the  sand  ejecta  using  the  present  model  with  parameters  for  Run6 

Turning  to  the  experimental  data  summarised  in  Table  2,  the  small  DOB  data  at  DOB  =  3  cm 
are  grouped  around  the  following  three  times,  101  psec,  126  psec  and  201-202  psec.  Therefore, 
plots  for  CTH  calculations  using  the  set-ups  of  Runs6-Runl0  are  only  shown  at  these  three 
times. 


Mesh  at  1  01c-04  seconds  Mesh  at  1  26c-04  seconds  Mesh  al  2  02e-04  seconds 


-20  -10  0  10  20  -20  -10  0  10  20  -20  -10  0  10  20 
X(cn»)  X(cm)  X(cm) 

Figure  19.  CTH  calculation  of  the  sand  ejecta  using  the  present  model  with  parameters  for  Run7 

As  in  the  previous  subsection,  we  start  CTH  analysis  using  the  present  implementation  with 
the  same  three  groups  of  the  yield  limit  parameters.  Numerical  results  for  these  parameters 
are  obtained  with  the  CTH  runs,  Run6,  Run7,  and  Run8.  Graphic  representations  of  these 
results  are  shown  in  Figs.  18, 19,  and  20  for  Run6,  Run7,  and  Run8,  respectively. 
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Mesh  al  1  Ole-04  seconds  Mesh  at  1  26c-04  seconds  Mesh  al  2  02e-04  seconds 
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X  (cm)  X  (cm)  X  (cm) 

Figure  20.  CTH  calculation  of  the  sand  ejecta  using  the  present  model  with  parameters  for  Run8 

Similarly  to  the  previous  calculation  series,  the  drastically  different  ejecta  shape  among  the 
results  obtained  with  the  CTH  calculations  employing  the  present  implementation  (Figs.  18- 
20)  is  observed  in  Fig.  18  for  the  case  of  the  high  yield  limits  at  both  compression  and  tension. 
It  is  interesting  to  note  that  the  apparent  ejecta  shape  appears  to  contradict  the  interpretation 
from  Runl  where  high  strength  at  tension  results  in  low  bubble  height.  However,  this 
apparent  mesh  contour  for  the  porous  material  should  be  looked  at  in  conjunction  with  the 
density  plots,  because  the  present  model  considers  a  two-phase  mixture  that  may,  for 
example,  contain  a  prevailing  volume  of  air. 


X  (cm)  X  (cm) 

Figure  21.  Density  distribution  inside  the  ejecta  for  Run6  calculation 


Indeed,  when  plotting  the  density  distribution  for  Run6  in  Fig.  21  along  with  the  material 
contours,  it  can  be  seen  that  the  high-density  (low-porosity)  material  is  concentrated  mainly  at 
the  lower  boundary  of  the  porous  material  and  directly  above  the  detonation  products 
pushing  the  ejecta.  Therefore,  the  actual  height  of  the  bubble  carrying  the  prevailing  mass  is 
much  smaller  and  it  is  adjacent  to  the  upper  boundary  of  the  detonation  products.  At  the 
same  time,  as  seen  from  Figs.  19-20,  the  prevailing  masses  of  the  ejecta  for  Run7  and  Run8  are 
concentrated  in  very  narrow  zones.  Thus  the  ejecta  contours  for  these  two  runs  approximately 
correspond  to  the  mesh  contour  of  the  ejecta  drawn  in  Figs.  19  and  20. 
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Mesh  at  1  Ole-04  seconds  Mesh  at  I  26c-04  seconds  Mesh  at  2  02e-04  seconds 
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Figure  22.  CTH  calculation  of  the  sand  ejecta  using  the  DRY  SAND  model  (Run9) 

The  final  two  CTH  calculations  of  the  present  series  were  conducted  with  the  DRY  SAND  and 
P-X  models  from  the  CTH  material  database.  The  numerical  results  of  these  two  runs  are 
summarised  in  Fig.  22  (Run9)  and  Fig.  23  (RunlO),  respectively. 


Mesh  at  1  01c-04  seconds  Mesh  at  I  26c-04  seconds  Mesh  af  2  02e-04  seconds 


X  (cm)  X  (cm)  X(cm) 

Figure  23.  CTH  calculation  of  the  sand  ejecta  using  the  P-X  model  (RunlO) 

Similarly  to  Run7  and  Run8,  the  calculation  with  the  DRY  SAND  model,  Run9,  has  the  ejecta 
mass  concentrated  in  a  very  narrow  zone.  Therefore,  the  ejecta  contours  in  this  case  also 
correspond  to  the  mesh  boundary  drawn  in  Fig.  22.  The  P-X  model  calculations  shown  in 
Fig.  23  result  in  an  uncertain  shape  of  ejecta  and  this  case  will  be  analysed  below  after 
comparison  with  the  experiment. 

Referring  to  the  experimental  data  listed  in  Table  2,  it  can  be  seen  that  out  of  the  3  tests  of  this 
series  (Shots  15, 16,  and  17),  the  results  of  Shots  16  and  17  are  more  consistent  to  each  other 
than  the  results  of  Shot  15.  Indeed,  the  bubble  height  for  snapshot  15-1  from  Table  1  is  close  to 
that  for  snapshot  17-1,  but  at  a  significantly  earlier  time  point.  While,  the  bubble  height  for 
snapshot  15-2  essentially  predominates  the  average  bubble  height  at  approximately  the  same 
time  points  from  Shots  16  and  17. 
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Figure  24.  Comparison  of  experimental  snapshots  at  late  stage  with  the  CTH  results 


Therefore,  for  comparison  of  the  numerical  results  at  a  late  stage  of  the  soil  expansion 
observed  in  the  experiments  we  take  the  x-ray  snapshot  16-2  that  is  quite  similar  to  the 
snapshot  17-2.  The  comparison  results  are  shown  in  Fig.  24  with  run  identifiers  for  the 
corresponding  CTH  calculation  in  each  frame.  From  the  comparison  it  might  appear  that  the 
bubble  height  of  Run6  agrees  better  with  the  experiment.  However,  as  commented  above  in 
Fig.  21,  the  two-phase  nature  of  the  modelling  should  be  remembered,  whereby,  the 
predominating  mass  bearing  contour  is  drawn  in  the  frame  corresponding  to  Run6  together 
with  the  external  ejecta  contour.  Comparison  of  this  part  of  the  ejecta  with  the  experimental 
contour  shows  that  the  constraining  effect  of  the  high  yield  limits  for  the  case  of  Run6  is  also 
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present  for  the  small  DOB  calculations.  Thus,  similarly  to  the  previous  series  of  calculations, 
the  last  two  runs  with  the  low  yield  limits  at  tension  (Run7  and  Run8)  produce  more  adequate 
results  when  compared  with  the  experimental  data. 

Another  candidate  for  a  good  description  of  the  bubble  height  is  the  P-X  model  calculation  of 
RunlO.  However,  when  referring  to  the  density  distribution  of  the  P-X  model  calculation 
shown  in  Fig.  25,  the  major  mass  of  the  ejecta  is  again  concentrated  in  the  lowest  part  of  the 
leading  piece  of  the  soil  ejecta  material  similarly  to  the  case  of  Run6.  The  corresponding  mass 
bearing  portion  of  the  calculated  ejecta  plotted  along  with  the  external  contour  of  the  ejecta  is 
shown  in  the  last  frame  of  the  comparison  graph  in  Fig.  24.  Thus,  the  contour  containing  the 
major  mass  of  the  ejecta  for  this  calculation  describes  the  experiment  in  Fig.  24  more  poorly 
than  the  similar  contour  for  the  case  of  Run6  criticised  earlier.  It  could  be  argued  that  the  fairly 
high  yield  limit  value  of  500  MPa  selected  for  this  calculation  constrains  the  bubble  height  of 
the  ejecta  in  the  present  case.  However,  when  referring  to  the  previous  subsection  calculations 
with  this  model  (frame  5  in  Fig.  16),  it  is  seen  that  a  reduction  of  the  yield  limit  would  result  in 
an  even  larger  discrepancy  for  that  case  shown  in  Fig.  16. 


Densitv  Mesh  at  2.01e-04  seconds 


x  (cm)  x  (cm) 

Figure  25.  Density  distribution  inside  the  ejecta  for  RunlO  calculation 

Thus,  from  the  comparison  of  the  experimental  snapshots  with  the  CTH  calculation  results 
shown  in  Fig.  24,  it  is  seen  that  Run8  results  are  reasonably  close  to  the  experimental  bubble 
shape  at  a  later  stage  of  the  soil  expansion.  When  correlating  the  present  model  calculation 
plots  with  the  yield  limit  input  data  from  Table  1,  it  can  be  seen  from  Fig.  24  that  a  larger 
ejecta  expansion  in  the  middle  section  of  the  bell-shaped  ejecta  for  small  DOB  is  likely  to  be 
associated  with  a  higher  strength  at  compression. 

Similarly  to  the  large  DOB  case  in  the  previous  subsection,  we  summarise  all  available 
experimental  snapshots  for  the  small  DOB  with  results  of  Run8  at  approximately  the  same 
time  points  in  Fig  26.  It  is  seen  that  the  CTH  calculation  using  the  present  implementation 
with  the  Run8  parameters  describes  the  experimental  data  (specifically.  Shots  16  and  17)  quite 
well.  As  commented  earlier,  possible  inconsistency  of  Shot  15  should  be  kept  in  mind. 
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Shot  17-2.  201.1  (is  3  cm  DoB  Shot  15-2.  201.9  (is  3  cm  DoB 

Figure  26.  Comparison  of  the  Run8  calculations  with  the  experimental  snapshots  [12] 


It  should  be  noted  that  even  good  agreement  of  the  ejecta  shape  with  the  experimentally 
observed  contours  does  not  guarantee  a  good  prediction  of  the  momentum  transfer  to  a  target. 
To  illustrate  this  point,  the  numerical  results  of  two  runs  for  the  small  DOB  at  t  =  50  psec  (one 
of  the  runs  with  the  present  model,  Run6,  and  Run9  with  the  DRY  SAND  model)  for  mesh, 
density,  pressure,  and  velocity  magnitude  are  superimposed  onto  4  single  plots  in  Fig.  27.  The 
left  half  of  each  graph  corresponds  to  the  calculation  with  the  DRY  SAND  model  and  the  right 
half  represents  the  results  calculated  with  the  present  implementation. 


UNCLASSIFIED 


DSTO-TR-2728 


Mesli  at  5.00e-05  seconds 


10 


X  (cm) 


Density  Density 


X  (cm) 


Pressure 


0 

X  (cm) 


Velocity 
I  10!  20 


Velocity 


0 

X  (cm) 


Figure  27.  Comparison  of  the  DRY  SAND  model  results  (the  left  half  of  the  plots)  with  the  present 
model  implementation  residts  (the  right  half  of  the  plots)  with  DOB  =  3  cm  at  t  =  50  psec 
(Run9  and  Run6  residts) 


It  can  be  seen  that  at  times  shown  in  the  graph  the  velocities  of  ejecta  are  approximately  the 
same  from  both  the  velocity  fringes  and  contour  positions.  The  contour  shapes  are  quite 
similar  at  this  moment  of  time  for  both  cases,  but  the  density  distributions  differ  dramatically 
(the  same  point  was  also  illustrated  by  the  one-dimensional  simulations  of  the  previous 
section).  The  sand  density  is  nearly  constant  through  the  ejecta  thickness  for  the  DRY  SAND 
model  with  the  SESAME  EOS,  whereas  the  material  compaction  in  the  present  model  is 
performed  in  a  kinetic  manner  demonstrating  a  significant  mass  gradient  through  the  ejecta 
thickness  (as  also  observed  in  the  one-dimensional  calculations).  Thus,  the  momentum 
transfer  can  be  significantly  overestimated  with  the  tabulated  SESAME  approach.  At  the  same 
time,  model  preference  cannot  be  established  before  verification  of  the  numerical  results  with 
experimentally  measured  through-thickness  density  distribution  or  with  recorded  momentum 
appropriately  resolved  in  time  has  been  conducted.  Appropriately  designed  experiments 
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should  allow  measurement  of  the  momentum  sensitivity  to  target  parameters  or  to  record 
directly  a  pulse  in  the  porous  sample  by  a  gauge. 


7.  Conclusions 


As  a  further  development  of  the  basic  CTH  implementation  capability  established  in  DSTO 
earlier  (see  [5]),  the  present  work  describes  another  CTH  implementation  of  a  two-phase 
porous  material  model  [4],  The  implemented  model  is  capable  of  describing  the  anomalous 
behaviour  of  highly  porous  materials,  keeping  the  formalism  consistent  with  that  of  the 
irreversible  extended  thermodynamics  (e.g.  [17]).  The  model  is  also  capable  of  describing  the 
phase  transition  effects  and  strain  rate  elasto-plastic  effects. 

The  present  implementation  is  consistent  with  the  CTH  user  implementation  interface  [3, 19]. 
This  means  that  1)  the  implementation  employs  the  mathematical  model  presented  as  the 
conservation  laws  with  a  set  of  constitutive  equations  evolving  along  the  particle  trajectory; 
and  2)  the  stress  response  is  decoupled  in  volumetric  (pressure)  and  deviatoric  (shear  stresses) 
responses.  The  implementation  is  realised  in  the  present  work  for  the  case  of  porous  materials 
and  the  implementation  flowchart  is  briefly  outlined  in  accordance  with  the  code  structure 
described  in  the  available  literature. 

The  model  implementation  has  been  verified  for  a  highly  porous  silica  material  against  the 
available  Hugoniot  data  [11]  and  against  calculations  using  an  in-house  code  [4,  6,  7],  The 
CTH  results  correlate  well  with  the  calculations  and  with  the  anomalous  behaviour  of  the 
material  observed  in  the  literature.  Driven  by  the  requirements  of  the  counter-mine  and 
counter-IED  project,  the  literature  data  [12]  on  the  explosion  of  a  100  g  C4  charge  buried  under 
a  layer  of  sand  have  been  simulated  with  the  present  model.  The  CTH  results  demonstrate  a 
good  description  of  the  experimental  data  by  accurately  following  the  soil  ejecta  shape. 

Results  of  the  present  report  have  demonstrated  that  an  advanced  CTH  implementation 
capability  for  multi-phase  models  has  been  established  in  DSTO.  However,  use  of  the 
implemented  models  concurrently  and  an  extensive  use  of  the  extra  variables  in  the  model 
implementation  will  require  even  more  intrusive  implementation  in  the  ERB  part  of  the  code 
along  with  the  possibility  of  routine  updating  of  extra  implemented  models  in  the  EOS  part  of 
the  code  (this  capability  is  currently  available  only  in  the  EP  part  of  the  code).  This  access  will 
also  be  critical  when  sharing  implemented  models  within  joint  projects. 

Addressing  the  accuracy  of  the  description  of  the  target  response  to  the  ejecta  effects,  further 
validation  is  necessary,  which  would  allow  one  to  evaluate  the  density  and  velocity 
distributions  through  the  ejecta  thickness  and  compare  them  with  experiments.  This  can  only 
be  done  by  validation  against  specially  designed  experiments. 
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